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Introduction 

This  Report  covers  the  results  of  modeling  work  performed  at  the  MIT  Space  Power  and 
Propulsion  Laboratory  for  the  past  2  1/2  years  on  two  Electric  Propulsion  technologies: 
Arcjets  and  Hall  Thrusters. 

The  Proposal  for  this  Grant  listed  two  Tasks: 

1.  Seeding  and  Multi-Species  Effects  in  Arcjets 

2.  Development  of  Two-Dimensional  Models  of  Hall  Thrusters 

This  was  supplemented  for  an  8-month  extension  with  proposed  work  on  Aicjet  Anode 
Attachment  and  on  Cross-Field  Transport  in  Hall  Thmsters. 

The  work  on  Task  1  (Seeding  and  Multispecies  Effects)  was  further  subdivided  into  studies 
of  Hj  segregation  in  N2-H2  arcjets  and  feasibility  studies  of  Cesium-seeded  hydrogen 
arcjets.  It  was  acknowledged  that  the  scope  of  proposed  work  exceeded  resources,  but  it 
was  hoped  that  additional  student  support  could  be  obtained  through  AASERT  program. 
In  the  event,  we  were  unable  to  secure  this  support,  and  the  subtask  on  H2  segregation  has 
not  been  performed. 

We  have,  on  the  other  hand,  performed  significant  work  on  Cesium-seeded  arcjets  and  on 
Hall  Thruster  two-dimensional  modeling.  This  work  will  be  briefly  described  here,  and 
further  documented  by  enclosed  Conference  papers.  The  work  on  the  Extension  areas  is 
stiU  incomplete,  and  will  be  separately  reported  at  the  end  of  a  planned  no-cost  extension. 

2.  Alkali-Seeded  Hydrogen  Arcjets 

The  motivation  for  this  work  is  the  existence  of  large  “frozen  losses”  in  conventional  arcjets 
using  molecular  gases.  These  losses  appear  to  be  unavoidable:  ionization  of  the  gas  is 

needed  for  conduction,  and  is  fairly  costly  in  terms  of  energy  (=15  eV  per  ion  for  usual 

fuels);  also,  at  the  temperatures  where  ionization  occurs,  molecular  dissociation  is 
complete,  further  adding  a  few  eV  per  eventual  ion.  Because  of  the  very  short  expansion 
time  and  the  moderate  pressures  in  use,  nozzle  recombination  is  very  limited,  and  most  of 
this  energy  is  lost  in  the  exhaust. 

Seeding  the  working  gas  with  small  amounts  of  an  alkali  (Cs  or  K)  has  been  standard 
practice  in  MHD  power  generation,  where  gas  temperatures  are  limited  below  ionizing 


levels  by  either  fuel  energy  (combustion  MHD)  or  wall  endurance  (non-equilibrium  MHD). 
Cesium,  in  particular,  has  an  ionization  potential  of  only  3.9  volts,  and  is  fully  ionized  at 
electron  temperatures  as  low  as  40(X)K.  This  could  be  directly  of  use  under  equilibrium 
conditions,  keeping  the  molecular  working  gas  neutral;  however,  significant  dissociation 
could  still  be  incurred  at  this  temperature.  A  second  physical  effect  can  be  exploited  at  this 
point:  since  the  Cs  fraction  will  be  kept  as  low  as  possible  to  avoid  overly  increasing  the 
gas  molecular  weight,  we  can  expect  a  relatively  low  conductivity  level,  and  hence,  for  a 
given  volumetric  heating  rate,  high  E  fields.  The  electron  temperature  is  known  to  increase 
with  E/P  (p=gas  pressure),  and  hence  we  can  expect  thermal  non-equilibrium  to  develop, 
with  Te  significantly  exceeding  the  gas  temperature  Tg.  Since  ionization  depends  on  Te, 
while  dissociation  is  mainly  a  function  of  Tg,  we  can  strive  for  conditions  where  Tg  is  no 
higher  than  perhaps  SOOOK  (sufficient  for  Isp  to  exceed  900  sec  in  hydrogen),  and 
dissociation  can  be  kept  to  a  minimum.  At  the  same  time,  Te  could  be  high  enough  to  fully 
ionize  the  Cs  seed,  although  still  below  some  level  of  the  order  of  7000K  which  would  be 
required  to  ionize  the  working  gas. 

A  secondary  benefit  of  the  increased  electric  field  level  is  higher  voltage  operation.  Within 
limits  (say  below  the  KV  level)  this  is  beneficial  in  that  it  minimizes  the  fractional  loss  due 
to  near-electrode  voltage  drops.  These  drops  are  themselves  likely  to  be  reduced  by  the 
presence  of  the  alkali,  perhaps  in  proportion  to  the  ionization  potential. 

One  final  effect  which  could  contribute  to  higher  efficiency  is  the  more  uniform  radial 
distribution  of  heating  in  a  seeded  gas.  The  strong  constriction  of  the  discharge  into  an  arc 
is  due  to  the  steep  increase  in  conductivity  associated  with  partial  ionization  on  the  arc 
edges;  in  a  situation  where  the  Cesium  can  be  fully  ionized  everywhere,  constriction  is  less 
likely.  This  has  two  consequences:  (a)  It  reduces  the  possibility  of  highly  localized 
regions  where  hydrogen  could  become  itself  ionized,  and  (b)  It  directly  improves  the 
efficiency,  because  in  any  propulsive  device  heat  is  most  efficiently  used  to  produce 
momentum  flux  if  distributed  uniformly  across  the  gas  stream. 

All  of  these  potential  benefits  must  be  eventually  weighed  against  (a)  Operational 
complexity  in  the  dispensing  of  Cs,  which  must  be  kept  liquid  and  its  vapor  must 
homogeneously  mix  with  the  hydrogen,  and  (b)  Concerns  about  plume  contamination. 
These  are  issues  for  later  study,  but  we  may  note  at  this  point  that  experience  with  Cs 
dispensers  already  exists  (ion  engines),  that  Cs  will  solely  be  used  in  about  1-3% 


proportion  by  mass,  and  that  Cs  has  a  high  enough  vapor  pressure  to  ensure  vacuum  self¬ 
cleaning  of  surfaces  at  temperatures  above  «270K. 

As  a  tool  to  assess  the  physical  parameters  and  the  performance  potential  of  a  Cs-seeded 
arcjet,  we  formulated  a  simple  constrictor  flow  model,  which  is  fully  described  in  paper 
mPC  95-236  (Ref.  1  appended).  The  model  contains  enough  2-D  information  to  allow 
computation  of  the  radial  profiles  of  Tg,  Te,  velocity,  etc.,  as  they  evolve  down  the 
constrictor.  It,  however,  cannot  account  for  near-electrode  effects,  and  it  also  ignores 
viscous  losses  (except  that,  to  keep  these  losses  reasonable,  the  constrictor  IVD  ratio  is  kept 
below  10).  The  discharge  is  assumed  to  attach  at  the  end  of  the  constrictor,  and  an 
isentropic  nozzle  expansion  is  then  assumed  to  occur  through  a  given  pressure  ratio.  This 
allows  a  simple  calculation  of  thrust  and  other  performance  parameters. 

The  paper  (Ref.  1)  shows  results  of  a  parametric  exploration  using  this  model.  The 
parameters  varied  were  the  seed  fraction,  chamber  pressure,  current  and  constrictor  inlet 
Mach  number  (this  last  parameter  controls  the  total  enthalpy  addition).  The  range  of 
variation  of  these  parameters  (one  at  a  time)  was  bounded  roughly  as  follows: 

— For  seed  fraction  below  =1.5%  by  mass,  the  electron  temperature  exceeded 

7500K,  because  of  the  high  fields,  as  noted  before.  This  was  taken  as  a  limit 
beyond  which  there  is  a  risk  of  ionization  of  the  hydrogen  itself,  potentially 
triggering  a  rapid  transition  to  a  conventional,  lower  voltage  arc  burning  in 
hydrogen. 

— For  seed  fraction  above  =2.5%  by  mass  the  length-to-diameter  ratio  of  the 
constrictor  exceeded  10,  probably  leading  to  high  viscous  losses. 

— For  chamber  pressures  below  2.8  atm.,  again  the  high  E/P  levels  led  to  Te 
being  over  7500  K. 

— For  pressure  above  4  atm,  L/D  exceeded  10. 

— For  currents  lower  than  16A,  L/D  also  exceeded  10 
— For  currents  higher  than  19A,  Te  exceeded  7500K 

— For  inlet  Mach  number  below  0.22,  Te  exceeded  7500K,  and,  simultaneously, 
the  dissociation  losses  increased  rapidly,  due  to  the  exit  gas  temperature  rise 


(This  was  manifest  in  the  corresponding  rise  of  Isp). 

— For  inlet  Mach  number  above  0.25,  Isp  became  lower  than  750  sec. 

All  these  results  apply  to  a  particular  choice  of  constrictor  radius  (1.4  mm).  The  central 
(nominal)  case  operated  at  10.1  KW  power,  V=595  Volt,  m  =  0.37g/x,  had  a  constrictor 
L/D  of  8.7,  and  delivered  a  specific  impulse  Isp=835  sec,  with  only  1%  frozen  losses.  A 
complete  calculation  of  efficiency  was  not  possible,  due  to  the  limitations  of  the  model,  but 
frozen  losses  generally  account  for  about  50%  of  the  input  power,  so  that,  even  allowing  a 
10%  viscous  loss  and  taking  no  credit  for  reduction  of  the  voltage  drops,  one  can  expect 
overall  efficiencies  approaching  80%. 

The  paper  also  presents  a  simplified  discussion  of  the  major  trends,  based  on  a  simple  heat 
addition  model,  made  possible  by  the  relatively  uniform  radial  distribution  of  power 
deposition. 

In  sum,  this  preliminary  exploration  has  yielded  very  encouraging  results,  and  we  plan  to 
pursue  the  idea  further.  Areas  where  obvious  advances  are  needed  are: 

(a)  Inclusion  of  viscous  effects 

(b)  True  2-D  modeling,  including  nozzle 

(c)  Simultaneous  modeling  of  hydrogen  ionization 

(d)  Stability  analysis,  particularly  against  ignition  of  a  second  (classical 
arc)  mode 

Other  areas  to  be  looked  into  are 

— Possible  segregation  of  Cs  away  from  center  of  discharge 

— Radiative  losses  (preliminary  estimates  show  small  losses, 
due  to  the  full  Cs  ionization) 

— Thruster  thermal  design 

— Design  of  Cs  dispenser 


— Plume  contamination 


3.  Two-Dimensional  Modeling  of  Hall  Thrusters 

Hall  thrusters  have  achieved  some  degree  of  technical  maturity,  and,  given  their  simplicity 
and  good  performance  in  the  middle  Isp  range,  are  currently  the  technology  of  choice  for 
many  missions.  At  the  same  time,  significant  improvements  are  still  possible  in  this  area: 
efficiency  could  be  significantly  raised,  flucmations  need  to  be  reduced,  plume  divergence 
remains  an  issue,  and  lifetime  extensions  would  be  welcome. 

In  contrast  with  the  extensive  empirical  basis  accumulated  in  Russia  over  the  past  20  years, 
there  is  still  a  lack  of  detailed  understanding  of  the  various  mechanisms  at  work  in  a  Hall 
thruster,  and,  especially,  a  lack  of  reliable  overall  models  which  incorporate  the  proper 
physics  and  use  it  to  project  performance  or  life.  This  has  been  our  motivation. 

An  earlier  effort  (Ref.  2),  using  a  1-D  model,  has  yielded  useful  performance  predictions 
capabilities,  although  it  lacked  any  ability  to  suggest  shape  improvements  or  to  predict 
erosion.  This  earlier  model  featured  particularly  the  assumption  that  the  electron  flow  to  the 
anode  was  throttled  by  the  radial  magnetic  field  following  Bohm’s  simple  cross-field 
diffusivity  formula.  Adequately  modeling  cross-field  diffusion  is  an  extremely  difficult 
task  in  plasma  physics,  and  the  success  of  the  1-D  model  using  the  simple  empirical  Bohm 
diffusivity  encouraged  us  to  develop  a  more  ambitious  2-D  version  still  using  this 
diffusivity. 

Collisionality  plays  a  very  minor  role  in  the  plasma  flow  in  Hall  thrusters  (in  fact,  the  onset 
of  significant  collisionality  limits  the  maximum  thrust  density  in  these  devices).  Thus,  the 
proper  theoretical  description  should  be  kinematic  rather  than  fluid.  On  the  other  hand,  full 
kinematic  modeling  of  the  very  fast  electron  dynamics  was  deemed  beyond  our  numerical 
capabilities,  except  possibly  for  very  simplified  geometries  (this  may  still  be  a  useful 
research  avenue,  as  shown  by  the  recent  work  of  Hirakawa  and  Yoshikawa,  Ref.  3).  It  is 
also  known  (and  this  was  early  called  Langmuir’s  paradox)  that  the  electron  population  in 
plasmas  is  far  closer  to  a  Maxwellian  distribution  that  would  seem  warranted  by  classical 
collision  considerations;  this  appears  to  be  due  to  the  randomizing  action  of  electrostatic 
fluctuations.  Because  of  this,  a  compromise  was  struck  in  our  model  by  adopting  a  fluid 
model  for  the  electrons  only,  while  using  a  Particle  in  Cell  (PIC)  description  for  ions  and 
neutrals. 

Whenever  fluid  descriptions  are  used,  the  price  to  pay  is  a  difficulty  in  defining  suitable 
models  for  transport  effects.  In  our  case,  in  particular,  cross-field  transport,  as  noted,  is 


the  main  issue,  and  it  is  fortunate  that  at  least  a  plausible  empirical  model  (Bohm’s)  is 
available.  As  will  be  shown,  this  ad-hoc  procedure  succeeds  remarkably  well  in  our  case; 
this,  however,  leaves  many  questions  still  unanswered,  such  as  the  role  of  wall  scattering, 
the  transition  to  classical  diffusivity,  etc.  These  will  be  commented  upon  later,  but  they 
remain  future  research  topics. 

The  model  and  its  capabilities  are  described  in  paper  lEPC  95-240  (Ref.  4,  appended).  It 
treats  an  arbitrary  axisymmetric  thruster,  using  a  precomputed  magnetic  field  from 
prescribed  polepiece  geometry.  Electrons  are  placed  in  Boltzmann  equilibrium  and  constant 
temperature  along  each  magnetic  line,  with  diffusion,  mobility  and  heat  conduction  across 
lines  of  force  which  are  derived  from  the  diffusivity  DBoh^=kTe/16eB.  Ions  and  neutrals 
are  lumped  into  50,000-100,000  macro-particles,  which  are  moved  at  each  time  step  in 
response  to  the  prevailing  electric  field  distribution.  This  field  is  updated,  together  with  the 
electron  properties,  between  ion-advance  steps.  Ionization  is  accounted  for  by  generating 
new  ion-electron  pairs  at  each  step  in  proportion  to  computed  ionization  rates.  Quasi¬ 
neutrality  is  assumed  everywhere,  so  that  electron  density  follows  directly  from  ion  density 
after  ion  advancement. 

The  model  was  exercised  for  the  conditions  of  the  only  existing  detailed  map  of  internal 
Hall  thruster  parameters,  taken  by  Bishaev  and  Kim  in  1978  (Ref.  5).  The  performance  of 
the  thruster  was  very  well  predicted  (Isp=1475  sec  vs.  1530  measured,  77  =  0.50  vs.  0.54 
measured,  I,„^g=3.10A  vs.  3.15  measured,  Ibe^=2.20A  vs.  2.10  measured). 
Significantly,  this  agreement  degraded  rapidly  when  the  Bohm  diffusivity  was  either 
increased  or  decreased  by  25%.  In  the  case  of  a  25%  decrease,  the  discharge  was  found  to 
essentially  quench.  It  is  to  be  noted  that  this  reduction  is  equivalent  to  a  similar  increase  in 
the  magnetic  field  B;  the  quenching  appears  to  be  the  result  of  a  reduced  ionization  rate  due 
to  too  few  neutralizer  electrons  penetrating  the  B  field,  and  of  the  subsequent  reduction  in 
current,  which  further  reinforces  the  trend. 

It  is  also  significant  that  the  model  was  able  to  calculate  accurately  the  ion  current  impinging 
on  the  insulating  walls  (1.1  Amp  vs.  1.0  measured).  Using  the  detailed  ion  impingement 
speeds  and  angles,  combined  with  experimental  sputtering  yields,  an  erosion  rate  of 

4  X 10^  m/ sec  was  calculated,  which  compares  reasonably  with  a  measured 
7  X 10^ m/ sec.  Finally,  the  “twist”  torque  due  to  side  magnetic  force  on  the  ions  was  also 


computed  and  found  to  be  3.1%  of  the  thrust  times  the  mean  radius.  No  data  exists  on  this 
quantity  for  the  Bishaev-Kim  conditions,  but  a  range  of  2%  to  8%  was  indicated  in  Ref.  6. 

The  detailed  comparison  of  internal  parameter  distributions  was  less  favorable,  although  the 
computed  results  were  qualitatively  correct.  In  particular,  the  peak  ion  density  was 
predicted  to  occur  closer  to  the  anode  than  observed,  and  to  reach  values  about  twice  too 
high.  This  may  be  associated  with  difficulties  in  modeling  correctly  the  neutral  injection 
pattern  form  the  complex  annular  anode-injector  used.  The  peak  electron  temperature 
(about  26  eV)  was  well  predicted,  but  the  data  showed  some  unmodelled  variation  along 
magnetic  lines,  which  may  be  due  to  magnetic  mirroring.  The  calculated  region  of  steep 
voltage  drop  (acceleration  region)  lies  mostly  outside  the  annular  enclosure,  while  the  data 
indicate  most  of  the  fall  to  occur  just  ahead  of  the  exit.  Perhaps  this  is  related  to  inaccurate 
representation  of  the  magnetic  field  shape,  for  which  insufficient  information  was  available. 

A  very  interesting  feature  of  the  computation  is  its  apparent  ability  to  mimic  the  observed 
unsteady  nature  of  the  flow.  In  fact,  all  of  the  properties  just  discussed  are  time  averages 
over  a  0.5  millisecond  window.  In  detail,  the  calculations  contain  fluctuations  of  current, 
temperature,  etc,  with  amplitudes  of  ±15%  and  with  a  variety  of  characteristic  frequencies. 
Of  these,  the  two  lowest  ones  (~6  Khz  and  28  Khz)  appear  to  be  related  to  neutral  transit 

time  through  the  channel  and  through  the  ionization  region,  respectively.  There  is  also  an 
80  Khz  component,  possibly  related  to  ion  transit  time,  and  components  above  100  Khz, 
which  may  be  contaminated  by  numerical  effects,  as  this  is  the  range  of  grid-crossing 
frequencies  for  the  particles.  Fluctuations  of  this  same  general  type  have  been  documented 
experimentally  in  many  SPT-100  tests.  Our  computations  appear  to  support  their 
interpretation,  advanced  by  V.M.  Gavryushin  (Ref.  7),  as  ionization  relaxation  effects. 

The  detailed  fluctuation  spectram  was  found  to  be  sensitive  to  details  of  gas  injection  and 
gas-wall  interaction  modeling,  which  is  reminiscent  of  reported  variability  with  operating 
conditions  and  with  thraster  aging. 

Following  these  computations,  efforts  have  been  made  to  better  characterize  the  critical 
cross-field  diffusivity.  Using  the  Bishaev-Kim  data  set,  the  electron  flux  was  backed  out, 
and  from  it  the  apparent  cross-field  diffusivity,  averaged  over  the  cross-section.  This 
appears  as  a  continuous  line  in  Fig.  1.  For  comparison,  the  Bohm  diffusivity,  using  our 
computed  B  field,  is  also  shown,  and  is  seen  to  be  too  high  over  the  downstream  half  of 


the  channel,  and  too  low  near  the  anode.  While  the  very  indirect  nature  of  the  data  make 
them  unreliable,  this  appears  to  indicate  approximately  the  correct  average  and  trend,  but  a 
defective  axial  distribution.  The  “1/B*”  diffusivity  curve  is  simply  a  curve  fit  used  for 
evaluation;  its  incorporation  into  the  model  did  yield  an  improved  position  of  the  peak 
ionization,  but  since  the  fit  lacks  any  physical  basis,  it  should  not  be  generally  adopted. 
Work  is  continuing  on  this  point;  there  are  indications  that  a  combination  of  Bohm  and 
classical  diffusivity  may  explain  the  rapid  rise  towards  the  anode  (where  electron-neutral 
scattering  is  strong). 

In  summary,  then,  our  Hall  thruster  hybrid  PIC  model  offers  exciting  insights  into  the 
physics  of  the  thruster,  as  well  as  a  clear  path  towards  applications  for  design.  Some  such 
applications  are  planned  for  the  near  future,  in  connection  with  microthruster  design  work 
under  different  sponsorship. 

At  the  same  time,  our  discussion  has  indicated  several  areas  of  model  deficiency  which 
need  to  be  addressed.  These  include: 

— Inclusion  of  electron-wall  scattering  as  a  cross-field  transport  mechanism 
— Inclusion,  similarly,  of  electron-neutral  scattering 
— Accounting  for  magnetic  mirror  effects 
— Improvement  of  the  magnetic  model 

— Improvement  of  the  gas-wall  model,  including  secondary  emission  details 

— Adaptation  to  the  TAL  geometry,  in  particular  the  very  short  ionization 
length  and  the  metallic  walls. 

Our  work  has  also  highlighted  the  need  for  additional  detailed  internal  data  in  Hall 
thrusters.  This  will  be  addressed  by  the  doctoral  student  responsible  for  this  work 
(Michael  J.  Fife)  under  a  separate  grant  from  NASA,  which  will  provide  access  to  its 
Lewis  R.C.  facility  and  to  an  operational  SPT-100  thruster. 


4.  Other  Work  Completed  during  the  Reporting  Period 

Although  initiated  under  an  earlier  Grant,  our  work  on  inlet  ionization  in  MPD  thrusters 
was  completed  during  this  reporting  period,  and  was  documented  in  Dr.  Eric  Sheppard’s 
Thesis  (Ref.  8)  as  well  as  in  papers  lEPC  93-218  and  AIAA  94-3342  (Refs  9,10 
appended).  This  work  showed  that  ionization  within  a  few  mm.  of  the  point  where  fresh 
gas  enters  the  MPD  channel  can  be  explained  on  the  basis  of  back-diffusion  of  charge  pairs 
from  the  dense  downstream  plasma.  It  also  established  conditions  under  which  convection 
will  overpower  diffusion  and  lead  to  discharge  “blow  off’. 

Also  completed  during  this  period  was  the  arcjet  model  development  by  Dr.  Scott  A. 
Miller.  This  was  documented  in  Dr.  Millers’  Thesis  (Ref.  10)  and  (partially)  in  paper  lEPC 
93-219  (Ref.  12,  appended).  The  bulk  of  the  work  had  been  reported  earlier,  in  our  1993 
Progress  Report,  as  well  as  in  Miller’s  paper  AIAA  93-2101.  This  work  advanced  in  a 
fundamental  way  the  state  of  the  art  in  arcjet  modeling,  by  accounting  for  electron  thermal 
non-equilibrium,  viscous  losses  and  a  variety  of  other  significant  effects.  Particularly 
important  was  the  finding  that  electron  temperature  elevation  explains  the  anodic  attachment 
of  the  arc  through  a  layer  of  cold  outer  gas,  with  no  need  for  artificial  conductivity  floors. 
In  Ref.  11,  Miller  also  coupled  the  gas  model  to  a  thermal  computation  of  the  arcjet  model, 
which  showed  only  moderate  feedback  on  performance,  but  allowed  rational  design  and 
life  estimations. 

Another  outcome  of  our  work  on  arcjet  modeling  was  a  review  by  this  author  of  this  whole 
field,  published  as  paper  AIAA  94-2653  (Ref.  13,  appended). 
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Figure  1:  Comparison  between  effective  electron  mobility,  Bohm  (1/B)  mobility,  and 
1/B®  mobility  in  the  acceleartion  zone  of  an  experiment^  SIT. 
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Abstract 

This  paper  traces  the  evolution  of 
arcjet  thruster  modeling  to  its  current  state 
of  development,  and  discusses  the  most 
pressing  theoretical  uncertainties  in  the 
way  of  further  progress.  The  thrust 
producing  mechanism  is  well  understood, 
and  so  is  the  general  energy  balance, 
leading  to  voltage  prediction,  provided 
anodic  attachment  is  separately  known. 
This  particular  problem  is  discussed  in  some 
detail.  It  is  argued  that  3-D  arc  root 
constriction  is  a  very  likely  cause  of  at  least 
some  of  the  observed  anode  behavior,  and 
should  be  theoretically  investigated. 
Other  issues  examined  are  electron- 
molecule  interactions,  species  segregation, 
and  numerical  efficiency  in  the  face  of  ever- 
increasing  complexity. 

1.  Introduction 

Like  many  other  practical  devices, 
arcjets  for  space  propulsion  have  evolved 
largely  on  the  basis  of  a  series  of  intuition- 
guided  empirical  developments,  with 
theory  providing  at  best  a  general 
background,  and  very  little  detailed  help 
to  the  designer.  While  scientific  data  and 
sub-models  have  existed  for  a  long  time 
regarding  the  properties  of  arcs  in  general 
(i.e.,  Ref.  1),  the  interactions  of  arcs  with 
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strong  flowfields  introduces  so  much 
additional  complexity  that  until  recently' 
there  was  little  hope  of  being  able  to 
construct  comprehensive  models  of  the 
complete  device.  The  situation  is  now 
rapidly  changing  due  to  the  combination  of 
greatly  expanded  computational  resources, 
and  renewed  interest  in  refining  arcjets  for  a 
variety  of  operational  space  uses.  In 
synergy  with  these  facts,  new  physical 
information  is  also  being  generated  via  a 
new  generation  of  plasma  diagnostics.  The 
net  result  is  an  exciting  rapid  increase  of  our 
understanding  of  arcjets.  Theory  still  lags 
substantially  behind  practice,  but  the  gap 
seems  now  bridgeable  for  the  first  time. 

In  this  paper  we  will  review  the 
arcjet  modeling  literature  (Sec.  2),  provide 
a  synoptic  summary  of  the  controlling 
physics  (Sec.  3),  examine  in  more  detail 
some  of  the  more  difficult  mechanisms  at 
work  (Sec.  4),  and  discuss  the  numerical 
methods  available  to  the  modeller  (Sec.  5). 
We  conclude  (Sec.  6)  with  a  look  ahead  and 
a  prognosis  for  the  near  future  in  this  field. 

2.  Literature  Review 

This  review  will  not  attempt  to  be 
comprehensive,  but  will  rather  be  used  to 
illustrate  the  main  lines  of  development  of 
this  area  in  the  past. 

2.1  Zero -Dimension  Models 

As  discussed  by  Pfender  in  his 
excellent  1978  review  on  arcs  and  arc 
heaters,  the  idea  of  using  an  electric  arc  for 
propulsion  purposes  grew  out  of  the 
previous  experience  with  arc  gas  heaters 
for  chemical  processing  and  other 
industrial  uses.  In  an  industrial  arc  heater 
the  object  is  a  uniformly  heated  gas 
exhaust,  and  devices  such  as  magnetic  arc 
rotators  have  been  used  to  homogenize  the 
temperature  in  a  gas  plenum.  Drawing  on 
this  analogy,  early  electrothermal  rocket 
models  evaluated  frozen  losses  (mainly 
dissociation  and  ionization  energy)  by 
postulating  uniform  conditions  in  the 
"heating  chamber",  from  which  the  gas 
was  assumed  to  then  flow  through  a  nozzle 
under  isentropic  (frozen,  or,  for  reference, 
chemically  equilibrated)  conditions.  When 
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applied  to  arcjets,  this  procedure  is 
fundamentally  flawed,  because  of  the 
combination  of  the  strong  segregation  of  hot 
and  cool  gas  in  an  arc  and  the  very  non¬ 
linear  variation  of  dissociation  and 
ionization  with  temperature.  Essentially, 
the  frozen  losses  are  confined  to  the  arc 
itself,  which  may  occupy  a  sizeable 
geometrical  fraction  of  the  arcjet  throat, 
but  carries  only  a  minor  part  of  the  exhaust 
flow  rate.  Nevertheless,  calculations 
based  on  such  models  did  call  attention  to 
the  large  losses  associated  with  frozen 
chemistry,  and  provided  at  least  a 
qualitative  explanation  for  the  loss  of 
efficiency  at  very  high  specific  impulses. 

2.2  Ouasi-One-Dimensional  Models 

After  the  initiation  of 
experimental  work  with  constricted  arcjets 
of  the  type  still  in  use,  the  essential  non¬ 
uniformity  of  the  flow  field  was  recognized 
and  accounted  for  in  models.  For  example, 
R.R.  John  et  al  subdivided  the  flow  into  a 
core,  a  transition  zone  and  an  outer  region, 
and  constructed  a  simple  analytical  model 
as  a  vehicle  for  interpreting  test  results  of 
their  30KW  arcjet.  This  work  identified 
the  basic  energy  addition  and  loss 
mechanisms,  including  friction,  chemical 
reaction,  area  ratio  effects  and  others. 
Similar  work  during  the  60's  and  70's  was 
reported,  for  instance  by  Weber  and  by 
Topham  1^1,  and,  more  recently,  by,  among 
others,  docker  et  al  and  Martinez- 

Sanchez  and  Sakamoto  and  Tahara  et 
al  [91 . 


The  persistence  of  this  type  of 
quasi-one-dimensional  modeling  despite 
the  increasing  availability  of  an  array  of 
more  elaborate  methods  is  due  to  their 
recognized  ability  to  mimic  the  most 
important  aspects  of  arcjet  operation  and  to 
provide  reasonable  performance  estimates. 
In  fact,  if  prediction  of  thrust  and,  to  some 
extent,  efficiency,  were  the  only  criteria  by 
which  models  were  to  be  chosen,  there 
would  be  very  little  incentive  for  the  vast 
increase  in  effort  required  at  the  next 
logical  modeling  level  (2-D  models). 


As  a  sample  of  the  quasi  1-D  type 
of  model,  a  brief  description  will  be  given 
here  of  that  of  Ref.  8.  The  flow  is 
subdivided  into  a  hot  core  of  radius  Ra  (the 
arc,  in  the  constrictor,  or  its  wake  in  the 
nozzle)  surrounded  by  isentropically 
expanding  uniform  gas,  out  to  the  radius 
R(x)  of  the  constrictor  or  nozzle.  The  arc 
grows  at  a  rate  determined  by  the  supply  of 
heat  radially  conducted  from  the  arc  core 

=  «> 

where  hae  is  the  "ate  edge"  enthalpy, 
corresponding  to  the  minimum  temperature 
for  electrical  conduction  (-70(X)K),  and  hout 
is  that  of  the  cold  outer  gas.  The  arc  carries 
a  total  flow  m^,  and  is  the  radial  heat 
flux  at  the  arc  edge.  The  static  pressure  and 
the  dynamic  head  q  =  pu^  are  assumed  to 
be  radially  uniform  (the  latter  based  on  its 
axial  growth  at  the  expense  of  p  decay,  if 
friction  is  ignored).  The  specific  enthalpy 
h  is  modeled  over  broad  ranges  of 

temperature  by  h  —  Q^p,  where 

Q=  y !  [y-\)  if  an  effective  ratio  of 
specific  heats,  y,  is  defined.  It  can  then  be 
shown  that  the  mass  and  enthalpy  fluxes 
are 


puh  =  -^jQpq^ 


(2) 


(3) 


and  the  total  enthalpy  flux  is  as  in  Eq.  (3), 
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with  an  extra  factor  .  -  ■  _  ^ 

2Qp. 

important  observation  from  (2)  is  that  most 
of  the  mass  flow  occurs  in  the  cool  layer. 
whereas,  from  (3),  most  of  the  power  flow 
occurs  in  the  hot  core.  This  has 
implications  for  overall  arcjet  behavior, 
which  will  be  discussed  later.  From 
relationships  such  as  (2)  and  (3)  above,  it 
follows  that  a  radial  profile  of  h(r/Ra)  is 
sufficient  to  specify  most  of  the  flow 
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variables.  This  profile  is  chosen  from  a 
solution  to  the  local  Ellenbaas-Heller 
equation  (radial  balance  of  Ohmic  heating 
and  heat  conduction).  The  thermal  and 
electrical  conductivities  are  correlated  to 
each  other  by  a  sort  of  Wiedemann-Franz 
law 


{^Kdr^B^C  (4) 

where  B  is  a  parameter  (2.79V  for  Hj, 
1.58V  for  N2).  After  this,  the  Ellenbaas- 
Heller  equation  gives  a  Bessel-function 
profile 
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and,  as  a  further  approximation, 
conductivity  and  enthalpy  are  also  related 
linearly  to  each  other: 


h^Gc  (6) 

"G  =  2.47xl0*/n'r2/sec"for  ^ 

2.20xl0^m^Q/sec^  foxNjj 

Given  these  functional  relationships,  it  is 
easy  to  construct  expressions  for  the  various 
integrated  fluxes  (mass,  energy,  current), 
and  to  derive  ordinary  differential 
equations  (ODE's)  expressing  their 
conservation.  The  resulting  set  of  ODE's 
has  a  singular  point  (sonic  passage)  which 
is  forced  to  be  at  the  constrictor  exit  by 
choice  of  stagnation  pressure,  given  the 
mass  flow  rate.  The  arc  attachment  point 
needs  to  be  externally  supplied,  and  is 
usually  placed  at  the  constrictor  exit.  The 
stagnation  temperature  of  the  outer  gas 
layer,  which,  as  noted,  carries  most  of  the 
flow,  is  an  important  parameter,  which 
needs  to  be  separately  estimated  by 
knowledge  of  the  thermal  configuration  of 
the  arcjet;  even  in  a  water<ooIed  case,  this 
estimation  is  made  difficult  by  the  heating 
effect  of  the  cathode.  By  contrast,  the 
starting  (cathode  end)  arc  enthalpy  is  an 
insensitive  parameter  which  seeks  its  own 
level  in  a  short  distance. 


The  degree  of  accuracy  obtainable 
from  this  level  of  modeling  can  be  judged 
from  Table  1,  where  comparison  is  made  to 
data  of  Ref.  10  (20Kw,  H2,  w'ater  cooled). 

As  noted  in  the  caption  in  Table  1, 
the  calculated  column  voltage  is  increased 
by  an  anode  drop  of  28V.,  after  which 
agreement  is  good.  In  hot-wall  cases,  the 
required  was,  however,  nearly  zero, 
which  may  not  be  realistic.  The  thrust  and 
related  quantities  (stagnation  pressure, 
specific  impulse)  are  well  predicted.  The 
"calculated"  heat  loss  is  simply  the 
assumed  anode  drop  times  the  current,  and 
is  shown  mainly  to  reinforce  the 
plausibility  of  the  value  chosen  for  AV^. 

Similar  levels  of  accuracy  are 
obtained  with  the  model  of  Refs.  6-7, 
which  also  require  choice  of  some 
adjustable  parameters. 

2A  Two-Dimensional  Models 

The  type  of  models  discussed  above 
have  serious  limitations  concerning  the 
details  of  the  gas  physics  and  flow 
dynamics.  While,  with  some  operator 
skill,  such  models  can  come  close  to 
providing  valid  performance  predictions, 
they  offer  little  hope  for  design 
improvement  based  on  exploitation  of  one  or 
other  of  the  many  competing  effects 
present.  The  inflexibility  arises  from  the 
need  to  postulate  simple  gas  physics  in 
order  to  construct  manageable  integrated 
fluxes,  whose  evolution  can  then  be  traced 
with  ODE's.  Relaxing  these  constraints 
requires  a  detailed  2-D  or  3-D  description 
of  the  problem. 

The  initial  work  on  2-D 
(axisymmetric)  arcjet  modeling  is  due  to 
Watson  and  Pegot  whose  focus  was  the 
constrictor  region  alone.  The  ga  “^assumed  to 
be  in  thermodynamic  equilibrium,  which  is 
reasonable  as  long  as  the  anodic 
attachment  region  and  the  nozzle  are  not 
examined.  Viscous  and  heat  conduction 
effects  were  included.  The  calculation  was 
made  by  spatial  marching,  using  what 
would  now  be  called  a  Parabolized  Navier 
Stokes  approach.  This  work  was  seminal 
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in  clearly  illustrating  the  interplay  of 
Ohmic  heating  and  flow  acceleration  and 
in  showing  in  detail  the  extreme 
nonuniformity  of  the  flow  field.  Efforts 
have  been  recently  made  to  supplement 
the  Watson-Pegot  model  with  near¬ 
electrode  and  nozzle  sections. 

Typical  of  the  state  of  the  art  about 
2-3  years  ago  are  the  models  of  Butler  and 
King,  of  Rocket  Research  Co.  and  of 
Rhodes  and  Keefer  1116]  ^f  UTSl.  The 
single-fluid  model  of  Butler-King  can 
handle  different  gases,  and  it  covers  the 
complete  thruster,  but  it  retains  the 
Watson-Pegot  equilibrium  assumptions.  As 
a  consequence,  the  electron  density  and  the 
electrical  conductivity  are  essentially  zero 
outside  the  arc,  and  no  simple  mechanism 
exists  for  current  to  reach  the  anode.  This  is 
remedied  by  introducing  a  "conductivity 
floor",  CT^,  of  the  order  of  0.1-1  (ohm)*l/cm. 
The  precise  O',  is  selected  to  provide  a  good 
prediction  of  total  voltage,  and  its  effect  on 
calculated  thrust  is  relatively  minor  (Fig. 
1,  from  Ref.  14).  This  is  analogous  to  the 
procedure  used  in  the  1-D  model  of  Ref.  (8), 
where  is  externally  supplied  to  match 
voltage,  with  minor  effects  on  thrust.  Also 
in  common,  once  the  constant  conductivity 
floor  is  imposed,  one  loses  the  ability  to 
predict  the  current  distribution  along  the 
anode  wall.  More  precisely,  the  model  does 
calculate  this  distribution,  but  there  is  no 
feedback  to  the  local  (7^  value;  as  we  will 
see,  this  feedback  leads  to  a  certain  amount 
of  axial  constriction,  which  is  here  ignored. 
These  difficulties  were  recognized  by  the 
authors  and  a  two-fluid  model  is  under 
development.  Cappelli  et  al  have 
compared  the  exit-plane  predictions  of  this 
code  to  data  for  a  1  KW  H2  arcjet.  Velocity 
is  well  predicted,  but  heavy-particle 
temperature  is  under-predicted 
significantly  on  the  jet  axis  at  the  higher 
powers. 


A  similar  situation  exists  with 
respect  to  the  model  of  Rhodes  and  Keefer 
^^^1.  This  model  incorporates  several 
physics  refinements,  such  as  a  calculation 


of  radiative  losses,  plus  a  full  set  of 
molecular  transport  effects.  However,  once 
again,  equilibrium  forces  a  very  low 
conductivity  outside  the  arc,  and  an 
artificially  elevated  value  is  used  for 
T<10,OOOK.  The  current  distribution  over 
the  anode  is  specified  a-priori,  which  for  a 
given  voltage,  results  in  excessive  Ohmic 
dissipation,  requiring  a  correction  factor  of 
0.75.  Of  course,  the  anode  is  not 
equipotential  in  general  when  current 
distribution  is  arbitrarily  prescribed. 


Variation  of  Pradicted  Current  with 
Electrical  Conductivity  Floor. 
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Variation  of  Predicted  Isp  with  Power  for 
Different  Values  of  Conductivity  Roor. 


Fig.  1:  Effects  of  conductivity  "floor”  (from 
Ref.  14) 


A  comparable  model,  with  a  single 
temperature  and  Saha  equilibrium,  has 
been  presented  by  Fujita  and  Arakawa  [17). 
Once  again,  an  ad-hoc  remedy  was 
necessary  in  order  to  generate  sufficient 
near-anode  conductivity,  and  it  took  the 
form  of  an  ionization  fraction  minimum  of 
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10'^.  The  single-temperature  assumption  is 
probably  even  less  adequate  in  this  case, 
because  Argon  was  the  test  gas,  and  heavy 
monoatomic  species  couple  very  p>oorly  to 
electrons  in  terms  of  energy. 

Other  single  fluid  models  have 
been  reported,  most  recently  those  of 
Andrennucci  et  al  and  Okamoto  et  al 
In  Ref.  (18),  the  focus  is  on  the  fluid 
dynamics,  and  Ohmic  heating  enters  as  a 
prescribed  input  term  in  the  energy  equation 
only.  In  Ref  (19)  the  temperature  and 
ionization  fraction  distributions  are 
peculiar,  vtith  maximum  values  ahead  of 
the  arc;  in  fact,  from  the  reported  current 
vector  maps,  no  arc  is  evident,  only  a 
convergence  of  current  lines  towards  the 
conical  cathode  tip. 

These  equilibrium  2-D  models  yield 
valuable  insights  on  phenomena  occurring 
in  various  parts  of  the  thruster,  but  care 
must  be  exercised  to  select  those  effects 
which  are  not  badly  affected  by  the 
artificial  correction  factors  introduced.  As 
far  as  overall  performance  prediction 
ability,  they  do  not  markedly  improve 
upon  the  simpler  1-D  models. 

The  first  model  to  abandon  the 
equilibrium  assumption  was  that  of  Miller 
His  model  endows  the  electrons  with  a 
separate  temperature  (Te^Tg,  thermal 
non-equilibrium  )  and  tracks  the  rates  of 
ionization  and  recombination,  rather  than 
using  Saha  equilibrium  at  Tg-  This  requires 
a  separate  electron  energy  equation,  which, 
aside  from  increasing  the  algebraic 
complexity,  also  requires  special 
integration  procedures  due  to  its  very  short 
characteristic  time.  Miller  also  included 
other  refinements,  such  as  ion-electron 
ambipolar  diffusion  and  finite  molecular 
dissociation  rates,  plus  a  very  detailed 
calculation  of  molecular  and  electronic 
transport  properties.  Flow  swirl  is 
accounted  for,  although  its  role  appears  to 
be  minor.  Only  Hydrogen  was  modeled, 
although  extension  to  pure  Nitrogen  is 
straightforward,  and  extension  to  N2-H2 
mixtures  is  possible,  but  laborious. 
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Fig.  2:  Radial  Profiles  of  Gas  and  Electron 
Temperature  0.25mm  Downstream  of  the 
Constrictor  Exit  for  the  Baseline  Case 
Arcjet  Simulation  (from  Ref.  31) 


0.00000  0.00040  0.00060  0.X120  0.00160  0.00200  0.00240  0.0O2SO 
R 


Fig.  3:  Radial  profiles  of  electron  density 
for  conditions  of  Ref.  20 
(m  =  O.lg  /  5,  /  =  60A),  at  several  axial 
sections  near  the  constrictor's  end. 

(a/  in  m~^,  R  in  m) 
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The  main  benefit  gathered  from 
these  improvements  was  the  ability  to 
calculate  the  detailed  current  path  outside 
the  arc  proper,  all  the  way  to  the  anode 
wall.  The  essential  point  is  the  strong 
thermal  nonequilibrium  existing  in  what 
we  could  call  the  "anodic  bridge"  region, 
namely,  the  annular  region  connecting  the 
arc  at  the  constrictor  exit  to  the  initial  part 
of  the  nozzle  wall.  This  is  illustrated  in 
Fig.  2,  which  belongs  to  a  simulation  of  the 
radiation  cooled  20KW  Hi  arcjet  of  docker 
etall^H-  In  the  arc  core  Te£T|despite 
the  very  large  Ohmic  power  input  to  the 
electrons;  this  is  because  of  the  high  degree 
of  ionization,  which  ensures  strong 
collisional  coupling  to  the  heavy  particles. 
But  when  we  move  outside  the  core  and  ng 
falls  rapidly,  this  coupling  is  weakened  to 
the  point  that  Tg  remains  in  the  10,000- 
20,000K  range,  while  Tg  asymptotes  to  the 
outside  gas  temperature  (-1000K).  This 
elevated  Te  makes  it  possible  for  the 
ionization  rate  to  keep  up  with 
recombination  and  diffusive  losses  and  to 
maintain  an  electron  density  ~10^*^m*^ 
right  up  to  the  wall  (Fig.  3).  This 

represents  an  ionization  fraction  of  only  10* 
but  is  sufficient  to  provide  the  required 
electrical  conductivity  in  this  critical 
region.  The  resulting  current  distribution  to 
the  anode  is  fairly  concentrated  in  the 
axial  direction  (Fig.  4).  One  word  of 
caution  needs  to  be  advanced  here:  Miller's 
model,  despite  its  success  in  self- 
consistently  determining  the  current  path 
downstream  of  the  constrictor,  is  still  forced 
to  model  the  constrictor  itself  as  an 
insulator:  otherwise,  the  current 

attachment  is  seen  to  creep  further  and 
further  forward  as  the  simulated  time 
progresses,  in  situatiorts  where  downstream 
attachment  is  known  experimentally  to 
occur.  We  will  return  in  Sec.  4  to  this 
problem,  and  will  also  use  Miller's  model  to 
illustrate  a  number  of  effects  and  point  to 
additional  unresolved  issues.  We  close  this 
discussion  with  some  performance 
prediction  results  (Table  2  and  Fig.  5)  from 
Miller’s  model,  again  compared  to  data  of 
Ref.  21.  The  degree  of  accuracy  is 
comparable  to  that  obtained  using  the  1-D 
models,  but,  of  course,  with  less  empiricism. 


Several  other  models  have  recently 
appeared  in  the  literature,  or  are  under 
current  development,  which  seek  to  build 
upon  the  lessons  so  far  learned,  in 
particular  the  need  for  non-equilibrium 
modeling.  This  includes  the  2-Fluid  version 
of  Butler  and  King's  model,  as  well  as  work 
by  Krier  et  al  and  Rhodes  and  Keefer 
Babu  and  Subramaniam  have 
initiated  work  on  detailed  accounting  for 
vibrational  nonequilibrium  in  molecular 
arcjet  working  gases.  This  effect,  on  which 
more  will  be  said  in  Sec.  4,  is  potentially 
important  in  N2-containing  gases,  such  as 
those  deriving  from  ammonia  on  hydrazine 
propellants. 

A  somewhat  separate  area  of 
modelling  is  that  concerned  with  plume 
flows.  Here  the  impetus  derives  from 
contamination  concerns,  as  well  as  from  the 
ready  availability  of  optical  and  other 
data  on  arcjet  plumes.  From  the  modeling 
viewpoint,  plume  flows  strongly  emphasize 
transport  and  non-equilibrium  effects,  even, 
as  Liebeskind  et  al  have  pointed  out  1^1, 
calling  into  serious  question  the  use  of 
continuum  models  (the  mfp  at  the  exit 
plane  of  their  IKW  arcjet  was  seen  to  be 
several  mm,  for  a  10mm.  diameter). 
However,  direct  velocity  measurements  at 
the  exit  plane  still  indicate 

insignificant  slip  between  species.  Boyd, 
Beattie  and  Cappelli  127,28)  have 
developed  and  validated  Monte-Carlo 
based  models  for  these  plumes  showing  a 
large  degree  of  rotational  temperature 
nonequilibrium.  In  a  related  effort, 
Fasoulas  et  al  l^^l  have  modeled  the  N2 
plume  of  a  hybrid  electrothermal- 
electromagnetic  thruster  used  at  Stuttgart 
Univ.  for  reentry  simulations.  In  this  case  a 
continuum  approach  was  used,  and  care  was 
taken  to  include  vibrational,  as  well  as 
electron  thermal  nonequilibrium.  The  N2 
low-density  chemistry  was  also  modelled, 
and  good  agreement  was  found  to 
experimental  data. 
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Fig.  4:  Current  Streamlines  for  the  Baseline 
Case  Arcjet  Simulation  (from  Ref.  31) 
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Table  2;  Comparison  of  Discharge  Voltages 
for  the  Cases  in  Fig.  5 


Fig.  5:  Comparison  of  Predicted  Sp>ecific 
Impulse  to  Experimental  Data  for  German 
TTl  Radiation-Cooled  Arcjet  (from  Ref.  31). 

3.  The  Essential  Ardet  Physics 

The  art  of  modeling  consists  of 
retaining  what  is  fundamental  for 
explaining  a  particular  effect,  while 
discarding  or  approximating  what  is  not. 
This  Section  will  attempt  to  provide  some 
of  this  persf>ective  for  arcjets.  Some  of  the 
points  to  be  made  are  well  established, 
others  are  little  more  than  informed  guesses 
based  on  as  yet  incomplete  ii\formation. 

3T.  JPigdiction  of  Thrust  and  Specific 
Impulse 

This  is  the  easiest  of  the  goals  of  a 
thruster  modeller.  In  the  case  of  an  arcjet, 
our  discussion  of  1-D  models  (Sec.  2.2)  can 
serve  as  our  starting  point.  It  was  noted 
there  that  very  little  mass  is  carried  inside 
the  arc;  in  essence,  the  arc  acts  as  a  fluid- 
dynamic  partial  plug  on  the  flow,  so  that, 
if  Aa  is  the  arc  cross  section  suitably 
defined  and  A  the  constrictor  cross  section, 
the  flow-carrying  area  is  A-Aa.  The  flow 
chokes  at  the  constrictor  exit,  where  this 
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area  is  least,  and,  if  ^  and  7,  are 
the  outer  flow  stagnation  pressure  and 
temperature, 

m- - : -  (7) 

^oul 

with  is  the  familiar  "characteristic 
speed"  of  a  rocket  chamber,  proportional  to 
being  the  gas  constant).  On 

the  other  hand,  the  also  familiar 
"mechanical"  interpretation  of  a  rocket's 
thrust  states  that  this  thrust  (F)  arises 
from  the  lack  of  a  reaction  surface  for 

OtOm 

in  the  area  of  the  throat  (times  a  nozzle 
coefficient,  cp,  of  order  1-2,  accounting 
principally  for  the  pressure  forces  on  the 
diverging  nozzle  bell): 

= ‘^F^o.ou.A  (8) 

Notice  that  the  whole  constrictor 
area  A  is  used,  because  P  is  constant  across 
it.  Aside  from  some  incidental  change  in  cp 
due  to  peculiarities  of  the  arcjet  flow,  Eq. 
(8)  is  as  it  would  be  in  the  absence  of  the 
arc,  but  Eq.  (7)  is  very  different.  The 
specific  impulse  Isp  is  given  by 

,  F  /,  \  A 

=  (9) 

m  'A-A^ 

a 

where  the  factor  in  parenthesis  would 
obtain  with  no  arc. 

Eq.  (9)  shows  the  crucial  elements 
for  successful  specific  impulse  prediction: 

(a)  The  stagnation  temperature  of 
the  outer  gas  needs  to  be  well  predicted.  In 
an  operational  sense,  this  highlights  the 
importance  of  regeneration  cooling,  or  of 
any  steps  taken  to  raise  the  temjjerature  of 
the  outer  gas  (this  may  not  necessarily 
imply  raising  the  anode  wall  temperature, 
because  the  thermal  boundary  layer  is  quite 
thin). 

(b)  The  size  of  the  arc  at  the 
constrictor  exit  needs  to  be  accurately 


predicted.  Since  the  arc  growth  is  a  result 
of  radial  diffusion  of  Ohmic  heat  into  the 
co-flowing  outer  gas,  the  implication  is 
that  the  thermal  diffusivity  and  the 
electrical  conductivity  must  be  modeled 
accurately.  This  may  not  require,  however, 
more  than  a  parametric  transport  property 
representation  which  is  faithful  "in  the 
large",  even  if  some  level  of  detail  needs  to 
be  sacrificed.  As  an  example,  the  arc  radius 
computed  by  the  crude  model  of  Ref.  8  is 
compared  in  Fig.  6  with  optical 
measurements  of  Ref.  10. 

Operationally,  the  factor  A/(A- 
Aa)  in  Eq.  (9)  implies  that  the  constrictor 
length  and  diameter  should  be  designed 
such  that  the  arc  comes  as  close  as 
thermally  allowable  to  the  constrictor 
wall.  Since  different  gases  and  different 
operating  points  lead  to  different  arc  radii, 
this  fact  needs  to  be  also  taken  into  account. 
For  example,  the  speed  of  sound  in  N2  is 
about  1/4  that  in  H2  at  the  same 
temperature,  giving  the  arc  4  times  longer 
growth  time  in  a  given  constrictor  length 
when  running  in  N2.  There  is  also  some 
difference  in  heat  diffusivity  (as 
characterized  for  example  by  B  in  Eq.  (4), 
but,  to  a  first  order,  a  satisfactory 
constrictor  for  H2  will  lead  to  premature 
attachment  if  run  in  N2,  as  noted  in  Ref. 
(21). 

Other  factors  affecting  thrust 
prediction  include  the  actual  velocity 
profile  in  the  arc  periphery,  and  the  effects 
of  wall  friction.  As  an  example  of  the 
latter,  consider  the  constrictor  H2  flow  in 
the  20  KW  arcjet  of  Ref.  (21).  The 
Reynolds  number  based  on  its  5mm  length  is 
-10^,  giving  a  total  laminar  boundary 
layer  thickness  of  about  0.17mm 
(corroborated  by  the  computations  of  Ref. 
(20))  and  a  net  viscous  drag  of  0.016N. 
There  is  additional  drag  in  the  nozzle,  but 

its  magnitude  per  unit  length  {npu^R 

decreases  rapidly,  so  that  the  net  viscous 
drag  is  of  the  order  of  0.05N,  compared  to 
thrust  levels  in  the  0.5-1 .5  N.  The  relative 
effect  will  be  somewhat  larger  at  smaller 
powers  (roughly  as  (Power)'^/^.  The 
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message  here  is  that,  while  important 
because  of  the  low  Reynolds  numbers 
involved,  arcjet  viscous  thrust  losses  remain 
secondary. 


Fig.  6:  The  experimental  radius  is  that  for 
1/e.  luminosity  decay  compared  to 
centerline(Ref.  8) 


As  interesting  as  what  is  essential 
for  thrust  prediction  is  what  is  not.  In 
particular,  Eqs.  (8)  and  (9)  show  that  the 
conditions  inside  the  arc  itself 
(dissociation,  ionization,  non-equilibrium, 
etc)  are  fairly  irrelevant  in  this  regards  as 
long  as  the  arc  radius  can  be  calculated. 
This  partly  accounts  for  the  success  of  many 
simple  models  in  thrust  and  specific 
impulse  predictiorts. 

3.2  Prediction  of  Voltage 

For  a  fixed  current,  voltage 
determines  the  power  absorbed  by  the  arc, 
and  its  calculation,  involving  as  if  does  an 
accounting  of  energy  losses,  is  generally 
more  difficult  than  that  of  the  thrust. 

The  voltage  is  naturally  divided 
into  the  column  voltage  and  the  electrode 
drops.  The  main  difficulty  with  column 


voltage  is  the  determination  of  column 
length,  i.e.,  the  location  of  the  anodic  arc 
attachment.  This  is  one  of  the  main 
unresolved  issues  in  the  field,  and  we  will 
return  to  it  in  Sec.  4.  If  the  length  is 
separately  known,  the  field  E=dV/dx  is  no 
more  difficult  to  calculate  than  the  arc 
radius  Ra-  This  follows  from  the  radial 
energy  balance 

(10) 

R. 


or,  using  O  » 


,  with  c  referring  to  the 


arc  centerline. 


.  JK(T,-T.)la, 

‘  K 


11) 


The  quantity  in  the  square  root  can 
be  recognized  at  the  of  Eq.  (4),  and,  in 
fact,  the  Ellenbaas-Heller  equation,  as 
solved  in  Ref.  (8),  yields  simply 


£  = 


2.45 


(12) 


The  phenomena  in  the  cathode  spot 
area  are  too  complex  for  the  level  of 
modeling  likely  to  be  of  interest  for 
propulsion,  but  at  the  same  time,  their 
performance  impact  is  sununarized  by  one 
insensitive  quantity,  namely  the  voltage 
drop  AVj,  for  which  empirical  information 
does  exits  [1],  [30).  A  separate  issue  relates 
to  cathode  lifetime,  which  has  so  far 
largely  defied  modeling,  because,  in 
addition  to  mass  loss,  one  would  also  need 
to  deal  with  morphology  changes  (whisker 
growth,  for  example),  which  tend  to  be  the 
life-limiting  factors  in  arcjet  applications. 

Accurate  prediction  of  the  anode 
fall,  AV^,  is  probably  bound  up  with  the 
issue  of  attachment  determination  (Sec.  4). 
On  the  other  hand,  the  magnitude  of  AV^ 
is  becoming  more  amenable  to  first- 
principles  calculations.  Two  components  of 
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AV^  can  be  discerned:  (a)  A  resistive  drop 
occurring  in  the  electron-deficient  (but  still 
quasi-neutral  and  collisional)  layer 
adjacent  to  the  anode.  This  layer  is  thin 
(-0.1mm)  and  is  determined  by  the  balance 
of  ambipolar  losses  to  the  wall  and  excess 
ionization  sustained  by  the  enhanced 
radial  field.  The  computations  of  Miller 
(20,31)  have  been  able  to  resolve  this 
sublayer  in  the  context  of  a  whole-thruster 
simulation.  Fig.  7,  from  Ref.  (20),  shows  a 
profile  of  voltage  from  the  cathode,  along 
the  arc  axis  and  across  the  anodic  bridge  to 
the  anode  surface.  The  rapid,  but  resolved, 
voltage  increase  near  the  anode  surface  can 
be  noted. 


where  ng  is  the  electron  density  at  the 
sheath  edge,  c,  is  the  mean  thermal 

velocity  of  electrons,  and  \g  =  ^[kT,  /  m- 
is  the  Bohm  velocity,  with  which  ions 
enter  the  ion-attracting  sheath.  Here  j  is 
the  local  current  density  at  the  anode 
surface,  which  needs  to  be  self-consistently 
computed  at  the  same  time. 

The  preceding  description  of  the 
anode  is  not  incompatible  with  the 
possibility  of  three-dimensional 
constriction  (see  Sec.  4),  but  actual 
resolution  in  a  3-D  computation  is  still 
beyond  reach. 


Fig.  7:  Electric  Potential  Axial  Profile  for 
the  Baseline  Case  Arcjet  Simulation  (Ref. 


31). 


(b)  A  sheath  voltage  drop, 
typically  negative  (anode  below  local 
plasma  potential),  which  is  required  to 
turn  back  the  excess  thermal  flux  of 
electrons  beyond  those  required  by  the 
circuit.  This  contribution  is  <0) , 

given  by 


j  = 


0.61ii£ 


(13) 


In  summary,  adequate  calculations 
of  voltage,  and  hence  efficiency,  are 
becoming  possible,  although  many  models 
are  still  incapable  of  direct  determination 
of  anode  drops.  This,  plus  uncertainties  in 
cathode  and  column  voltages,  typically 
yields  scatter  of  +20-40V  in  reported 
results. 

3.3  Prediction  of  Heat  Losses 

This  is  important  from  the  thruster 
integration  point  of  view,  in  addition  to  its 
direct  effect  on  efficiency.  The  dominant 
contributor  to  wall  heat  flux  is  simply  the 
product  of  anode  drop  times  current,  with 
ordinary  conduction/convection  playing  a 
secondary  role.  While  this  somewhat 
facilitates  the  p>erformance  modeling  task, 
it  also  means  that  an  accurate  prediction  of 
the  heat  flux  distribution,  and  hence  of  the 
peak  wall  temperature,  depends  critically 
on  being  able  to  calculate  the  current 
distribution  at  the  anode,  which,  as  we 
have  seen,  is  still  somewhat  elusive.  Once 
the  wall  flux  is  determined,  it  is  a 
relatively  easy  task  to  calculate  the 
temperatures  throughout  the  anode  block, 
accounting  for  either  radiation  or 
regenerative  cooling.  One  caveat  in  this 
respect  is  that  the  wall  thermal  reaction 
time  (  a  few  seconds)  is  orders  of 
magnitude  longer  than  the  flow  reaction 
time  (tens  of  microseconds),  so  that  a  fully 
coupled  computation  is  impractical. 
Fortunately,  as  has  been  shown  by  Miller 
(31),  the  detailed  wall  temperature 
distribution  has  only  a  minor  effect  on  both. 
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performance,  and  current  distribution  in  the 
plasma  (as  long  as  the  gas  temperature  is 
nearly  correct).  It  appears  that  one  or  two 
cycles  of  iteration  on  the  wall  temperature, 
starting  from  a  reasonable  guess,  are 
sufficient  for  convergence. 

As  an  example,  the  heat  flux 
distribution  for  three  conditions  from  Ref. 
(31]  are  shown  in  Fig.  8.  Here  the  wall 
temperature  was  pre-assigned,  of  the  order 
of  HOOK.  With  these  fluxes  as  inputs,  a 
computation  of  the  heat  flow  in  the 
thruster  body,  including  radiation  boundary 
conditions  on  the  outside,  and  partial 
regeneration  to  the  incoming  hydrogen,  as 
in  the  experiments  of  Ref.  (21),  yielded  the 
wall  temperatures  of  Fig.  9,  with  a  hot  spot 
in  the  anode  attachment  area.  Repeating 
now  the  flow  computation,  with  account  for 
these  wall  temperatures,  yielded  almost 
the  same  heat  fluxes,  and  modified  the 
sp>ecific  impulse  to  the  extent  shown  in  Fig. 
10. 


Fig.  9  Anode  Wall  Temperature  at 
m  =  0Ag/s  for  Three  Applied  Currents 
(Ref.  31) 


_ Is60A 

_ IslOOA 

. I=130A 


Fig.  8:  Heat  Flux  to  anode  at  m  —  OAg/s 
for  three  applied  currents  (Ref.  31) 


Fig.  10:  Effect  of  Anode  Model  Coupling  on 
Predicted  Performance  for  m  =  Q.lg  / S 
(Ref.  31) 
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4.  f^ome  Modeling  Challenges 

The  preceding  section  may  seem  to 
imply  complacency  about  the  current  status 
of  arcjet  modeling.  But,  in  fact  if  we  did  not 
have  the  empirical  data  base  we  have 
about  the  current  standard  arcjet  design,  it 
is  doubtful  we  would  be  able  to  "invent  it" 
in  software.  Similarly,  with  the  current 
modeling  capabilities  it  is  not  easy  to  break 
away  from  the  beaten  path  and 
realistically  evaluate  new  designs 
(including  such  apparently  minor 
variations  as  alkali  seeding  or  subsonic 
attachment).  The  reason  is  that  there 
remain  several  basic  effects  which  are 
poorly  understood,  and  others  where  the 
basic  scientific  data  exist,  but  have  yet  to 
be  incorporated  in  our  models.  This  Section 
will  briefly  review  the  most  salient  of 
these  issues. 

4..1  Arc  Attachment 

As  implied  in  the  discussion  so  far, 
arc  attachment  at  the  anode  looms  large  as 
a  problem  area.  To  different  degrees,  all 
existing  models  are  forced  to  incorporate 
external  information  about  it  in  order  to 
proceed. 

In  the  physically  detailed  model 
of  Millert^’^^'l^^ '  there  remains  a 
troublesome  tendency  for  the  arc 
attachment  to  migrate  upstream,  which  is 
checked  by  modelling  the  constrictor  itself 
as  an  insulator.  It  will  be  useful  to  examine 
more  closely  the  assumptions  used  in  that 
model,  as  a  basic  for  fUjjher  discussion.  As 
noted  in  Sec.  2,  ionization  is  kinetically 
controlled  using  the  equation  (in  steady 
state). 

=  (14) 

where  includes  an  ionization  rate 
(strongly  dependent  on  Te)  and  a 
recombination  rate  (depending  on  ne).  The 
ion  velocity  V-  is  the  sum  of  the  mean  mass 
velocity,  u,  plus  the  "diffusion  velocity" 
V,,  and  in  this  model,  n,V,  is  w'ritten  as 


—Dyn,,  where  Da  is  the  ambipolar 
diffusivity.  Furthermore,  it  is  assumed 

that  only  the  radial  component  of  V,.  is 
strong  enough  to  warrant  consideration. 
The  electron  temperature,  Te,  which 
controls  ionization,  is  found  from  a  detailed 
local  energy  balance,  including  Ohmic 
heating,  ionization  work,  collisional 
cooling  by  heavy  particles  and  electron 
heat  conduction.  As  shown  before,  the 
heavy  gas  temperature  (controlled  by  its 
own  energy  equation)  and  the  gas  velocity, 
are  found  to  be  nearly  undisturbed  in  the 
anode  region  despite  the  current  bridge. 

The  picture  that  emerges  from  this 
set  of  assumptions  and  partial  results  is 
that,  in  the  2-D  steady  state,  neutral  gas 
enters  the  anodic  bridge  at  about  sonic 
velocity,  and  is  subjected  to  the  ionizing 
action  of  the  overheated  electron 
population  crossing  this  bridge.  The 
position  of  the  "leading  edge"  of  the 
resulting  ionization  region  results  from  the 
combination  of  forward  electron  heat 
diffusion,  and  the  finite  ionization  delay 
time  implicit  in  Eq.  (14).  This  last  effect  is 
probably  not  controlling,  because  a  slight 
additional  Te  elevation  can  rapidly 
multiply  the  ionization  rate  to  the  level 
required.  The  "trailing  edge"  of  the  anodic 
bridge  area  results  similarly  from 
downstream  electron  heat  diffusion  and  net 
recombination  once  Te  falls  below  some 
minimum  value.  Here  the  kinetic  delay 
does  play  a  role,  because  Te  does  not  control 
recombination  rate  and,  there  is  a  diffuse 
recombination  trail,  in  contrast  to  a  fairly 
abrupt  ionization  leading  edge.  The  point 
to  be  made  about  this  leading  edge  is  that 
the  electron  energy  equation,  which 
determines  the  diffusive  spread  of  the 
high  Te  region,  does  not  contain  a  strong 
downstream  convective  term,  as  the  Tg 
equation  does;  wherever  the  dissipative 
bridge  may  be  located  at  a  given  time,  a 
certain  region  ahead  of  it  receives 
electronic  heat  conduction  —K^T^,  and 
the  bridge  tends  to  move  forward. 
Incidentally,  the  same  cannot  be  stated 
about  the  downstream  "edge",  because  the 
dissipation  rate,  which  feeds  Te,  decays 
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rapidly  once  the  flow  divergence  spreads 
the  arc. 

It  appears,  therefore,  that  in  an 
anodic  bridge  sustained  only  by  hot 
electrons,  with  no  heavy  gas  disturbance, 
there  is  no  obvious  effect  capable  of  pushing 
the  bridge  downstream.  Either  a  surface 
effect  is  required  that  will  anchor  the  arc 
foot,  or  some  stronger  involvement  of  the 
heavy  gas  must  occur.  The  first  of  these 
alternatives  is  unlikely,  because,  unlike 
the  cathode,  the  anode  is  a  fairly  passive 
structure  in  general.  As  to  the  second 
alternative,  if  the  heavy  gas  temperature 
remained  fairly  well  coupled  to  the  gas 
temperature,  convective  transport  would 
then  imply  counterflow  heat  diffusion  inthe 
frontal  part  of  the  anodic  bridge,  with  a 
clear  trend  to  sweep  the  bridge 
downstream.  Let  us  examine  this  point 
more  closely. 

The  relative  strength  of  the 
various  terms  in  the  electron  energy  balance 
in  the  anode  attachment  region  for  the 
baseline  case  (rfi  =  O.lg /sec,  /  =  lOOA) 
of  Refs.  (20,31)  is  illustrated  in  Fig.  11,  from 
Ref.  (31).  The  input  is  clearly  Ohmic 
heating,  and  the  energy  sinks  are 
ionization,  dissociation  and  collisional 
loss,  with  similar  contributions  (although 
ionization  dominates  very  near  the  wall). 
It  would  therefore  appear  that  we  have  a 
borderline  situation,  such  that  an  increase 
in  density,  which  favors  collisional 
transfer,  should  bring  about  a  significant 
heavy  temp>erature  increase  in  this  region. 
This  was  investigated  by  considering  a 
range  of  m  from  0.05  to  0.15  g/sec,  at 
constant  (Power/ rfj ),  such  that  the  main 
effect  is  a  density  change.  What  was  found 
(Fig.  12)  is  that,  as  density  increases,  the 
arc  constricts  more,  in  response  to  a 
decreased  heat  diffusivity  and  the  whole 
attachment  region  spreads  further 
downstream  to  lower  density  regions.  The 
gas  temperature  remains  nearly 
undisturbed,  and  the  electron  temperature 
continues  to  reach  the  20,000- 30,000K  level 
near  the  anode. 


Fig.  11:  Radial  Profiles  of  Some  Terms  in 
the  Electron  Energy  Equation  in  the  Current 
Attachment  Region  Just  Beyond  the 
Constrictor  Exit  (Ref.  31) 


Fig.  12:  Enclosed  Current  Contours  for 

—  ^WSMJIkg:  /  =  60A,  m  =  0.05^/5 

m 

(top); 

/  =  lOOA,  mO.lOg /5  (middle);  and 
7  =  130 A,  m  =  0.\5gls  (bottom)  (Ref.  31) 
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Perhaps  one  should  also  investigate 
variations  in  the  collisional  coupling  model 
(which  is  based  on  the  <5  factors  quoted  by 
Sutton  and  Sherman  as  a  possible  route 
towards  higher  Tg  and  more  convective 
effects. 

4.2  AzimuthalJv  Non-Piffuse  Attachment 

Closely  allied  with  this  subject  is 
the  issue  of  whether  the  axial  symmetry 
assumption  is  valid.  It  is  conceivable  that 
the  arc  bridge  might  become  three- 
dimensionally  constricted  at  some  point, 
resulting  in  attachment  via  "spokes" 
(which  would  rotate  due  to  the  flow  swirl). 
In  that  case,  the  much  higher  local  plasma 
density  in  the  spoke  would  indeed  elevate 
Tg  (as  in  the  arc  core  itself,  see  Fig.  2)  to 
near  Te,  and  the  attachment  physics  would 
radically  change.  In  particular,  models 
such  as  those  of  Maecker  and  its  follow- 
on  by  Yamada  et  al  now  become 
relevant.  Maecker  argues  that  an  arc  moves 
by  a  superposition  of  mass  motion  of  the 
heated  gas  and  motion  of  the  high 
temperature  zone  with  respect  to  the  gas. 
In  a  steady  arc  in  cross  flow,  such  as  the 
anodic  bridge,  these  contributions  add  up  to 
zero  velocity.  Maecker  proceeds  to  relate 
the  local  curvature  of  the  arc  to  the  cross- 
flow  velocity  (as  in  a  horizontal  arc  subject 
to  upwards  convective  flow),  and  Yamada 
et  al  extend  the  model  to  calculate  the 
point  of  attachment  to  a  cylindrical  anode. 
Aside  from  the  relative  crudity  of  some  of 
the  assumptions  made  (for  example,  the 
cross-flow  speed  in  the  arc  core  is  set  equal 
to  the  external  normal  speed,  times  the 
inverse  density  ratio,  as  if  the  flow 
direction  remained  undisturbed),  one  needs 
to  keep  in  mind  that  the  basic  Maecker's 
model  is  based  on  a  "bent  core",  with  high 
gas  temperature.  Three-dimensional 
constriction  appears  necessary  for  this. 

This  problem  of  the  stability  of  the 
axisymmetric  solution  remains  largely 
unexplored.  Miller  performed  some 
linearized  electrothermal  stability 
analysis,  assuming  only  planar  disturbances 
with  an  axial  wave  vector.  Although  the 
axisymmetric  stability  problem  would 
involve  tangential  wavevectors,  the 


difference  may  be  small,  because,  as  noted, 
convective  effects  play  a  minor  role  in  the 
basic  solution.  The  results  indicated 
stability  for  the  conditions  of  the  numerical 
studies,  but  they  also  showed  that 
instability  would  result  at  higher  current 
densities  and  lower  ionization  fractions; 
given  the  inaccuracies  in  the  stability 
model  (uniform  background,  undisturbed 
heavy  species),  the  potential  for 
constriction  cannot  be  dismissed.  Fig.  13,  for 
example,  maps  the  stability  region  as  a 
function  of 
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Fig.  13:  Neutral  Stability  Lines  as  a 
Function  of  the  Zeroth  Order  Electron 
Temperature  and  Current  Density  with 

Nonequilibrium  Ionization  (Ref.  31) _ 

normal  current  density.  Modes  with  wave 

number  <  5000m~'  (A  >  1mm)  can 
become  unstable  when  Te-25,000K  if  the 
basic  currents  density  exceeds  lO^A/m^;  the 
approximate  peak  current  density  in  the 
baseline  case  (0.1  g/sec,  lOOA)  was  3x10^ 
A/rrA. 


More  refined  stability  studies  of 
the  attachment  region  are  desirable,  given 
the  important  implications  of  constriction. 
These  studies  are  analytically  challenging, 
and,  with  steady  progress  in  the  numerical 
area,  perhaps  a  direct  3-D  numerical 
attack  is  a  near-term  possibility.  An 
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interesting  intermediate  course  would 
involve  numerical  solution  of  a  linearized 
3-D  model,  starting  from  a  converged  2-D 
simulation,  although  any  linear  analysis 
will  miss  the  constricted  end  state,  which 
can  only  be  captured  by  a  full  3-D 
simulation. 

Some  experimental  evidence  of  3-D 
constriction  in  the  attachment  area  has 
been  reported  by  Harris  et  al  1^51,  who 
tested  a  radially  segment  arcjet  and 
observed  generally  non-symmetric 
distribution  of  current  to  the  segments,  with 
random  fluctuations,  indicative  of  a  multi¬ 
filament  structure  of  the  anodic  bridge. 

Similar  experiments,  but  with 
axial  segmentation  of  the  attachment 
region,  were  performed  by  Bems  et  al 
In  this  case,  the  attachment  was  confined  to 
an  extended  cylindrical  region  upstream  of 
the  nozzle  throat,  as  in  the  older  Giannini 
arcjet  design  with  arc  stabilization 
being  provided  mainly  by  the  flow  swirl. 
Berns  et  al  observed  oscillations  in  the 
segment  currents  and  overall  voltage  which 
indicated  downstream  convection  of  the 
attachment  point  until  a  new  re-strike 
occurred  near  the  cathode,  presumably 
when  the  local  arc-to-anode  voltage  in 
that  region  exceeded  some  threshold.  This 
same  re-strike  behavior  had  also  been 
mentioned  by  Pfender  ^  as  being  common  in 
plasma  troches.  In  Bern  s'  study,  when  the 
arc  foot  passed  through  the  throat  to  the 
supersonic  region  of  the  nozzle,  the 
operation  tended  to  be  re-strike  free,  but 
the  incidence  or  absence  of  re-strike  is  not 
well  documented  for  the  usual  supersonic 
attachment  designs.  In  any  case,  their 
occurrence  in  the  subsonic  attachment  case 
appears  to  lend  support  to  the  view  that 
the  anodic  bridge  may  be  constricted  in  a  3- 
D  sense,  because,  as  noted,  this  should  then 
lead  to  its  convection  by  the  flow. 

4.3  Molecular  and  Molecule-Induced 
Nonequilibrium 

A  very  different  area  where 
modeling  needs  improvement  is  that  of 
nonequilibrium  molecular  effects.  These 
include  slow  kinetics  or  freezing  of 


vibrational  and  rotational  modes  in  the 
nozzle  expansion  (certainly  in  the  plume), 
as  well  as  distortion  of  the  electron  energy 
distribution  associated  with  strong  cross- 
sections  for  inelastic  electron-molecule 
collisions. 

Of  these  two  issues,  the  first  is 
probably  the  least  significant,  although 
uncertainties  remain.  Zube  and  Myers 
were  able  to  measure  electronic  excitation 
temperature,  as  well  as  rotational  and 
vibrational  temperature,  in  a  IKW  H2*N2 
arcjet,  at  several  points  in  the  supersonic 
nozzle.  They  found  that  the  excitation 
temperatures  remained  above  10,000K  at 
the  nozzle  exit,  but  the  apparent  indication 
of  a  very  high  exit-plane  Te  was  dispelled 
by  a  more  refined  analysis  by  Zube  and 
Auweter-Kurtz  using  a  method  by  Park 
(40)  to  correct  for  lack  of  Boltzmann 
equilibrium  among  electronic  states.  The 
true  exit  plane  Te  was  found  to  be  3000- 
4000K,  in  agreement  with  model  results  by 
Miller  (®^(  and  with  measurements  by 
Hoskins  et  al  (41 J.  More  important  for  the 
present  discussion,  the  vibrational 
temperature  was  found  in  Ref.  (38)  to  track 
closely  the  rotational  temperature  during 
the  nozzle  expansion,  both  of  them 
decaying  to  2500-3000K  at  the  exit  plane. 
In  contrast,  in  a  5-1 OKW  N2  arcjet  study, 
Tahara  et  al  observed  that  the 
rotational  temperature,  after  an  initial 
drop  from  an  already  high  9000K  to  6000K, 
subsequently  rose  in  the  nozzle  to  an  exit 
plane  value  of  7OOO-10,O0OK.  These  tests 
were  at  fairly  low  pressure  (a  6mm 
diameter  constrictor  was  used,  giving 
chamber  pressures  of  0.1  to  0.3  atm),  which 
may  explain  the  anomalous  behavior. 

In  terms  of  modeling,  accounting  for 
vibrational  lag  is  relatively 
straightforward,  as  in  Ref.  (29)  for 
example,  although  it  does  involve 
additional  computatiorwl  burden.  A  more 
comprehensive  approach,  in  which  the 
stale-specific  vibrational  populations  are 
tracked  separately  (2^<<(^2!^  by  contrast, 
quite  demanding  in  the  context  of  a  2-D 
computation,  and,  in  view  of  the  data, 
would  not  seem  justified  for  the  normal 
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operational  regimes.  It  may,  however,  be 
needed  in  order  to  understand  the  effects  on 
the  free  electron  distribution,  to  which  we 
now  turn. 

It  is  well  known  that  molecular 
gases,  particularly  N2,  can  readily  absorb 
collisional  energy  from  the  free  electrons 
into  their  vibrational  modes.  For  N2,  this 
occurs  mainly  for  electron  energies  from  2  to 
3  eV.  As  a  result,  for  values  of  the 
field/density  (E/N)  ratio  high  enough 
that  significant  population  of  electrons 
would  exist  above  these  energy  levels  in 
atomic  gases,  very  few  actually  exist  in  a 
molecular  gas.  The  distribution  "dams  up" 
agaiirst  the  effective  barrier  represented  by 
the  vibrationally  active  region  of  the 
energy  spectrum,  and  it  greatly  departs 
from  the  Maxwellian  shape.  Typical 
results  by  Nighan  are  showm  in  Fig.  14. 
For  reference.  Fig.  15  shows  calculated  E/N 
in  the  anodic  bridge  region  of  a  H2  arcjet, 
using  the  model  of  Refs.  (20,31). 

For  N2,  if  £ /N  <  5  xlO’^Kcm^, 
there  develops  a  substantial  depletion  of 
high  energy  electrons.  If  this  persists  all 
the  way  to  the  ionizing  range  (13'16  eV), 
the  ionization  rate  calculation  would  be 
substantially  affected.  In  this  connection, 
it  must  be  noted  that  for  each  E/N,  the 
results  of  Nighan  appear  to  indicate  that 
the  distribution  actually 


Fig.  14;  Electron  energy  distribution 
functions  in  N2  for  various  E/N  values.  The 
distribution  function  is  defined  such  that 


J  u^f{u)du  =  1  and  the  reduced  average 

2  r  y 

energy  such  that  u'^f{u)du. 

(Ref.  43) 


Fig.  15:  Radial  Profiles  of  E/N  (V  cm^)  for 
conditions  of  Ref.  21 

(0.1  g/s,  60  A)  at  several  axial  sections  near 
the  constrictor’s  end. 


becomes  super-Maxwellian  at  high  enough 
electron  energy;  however,  the  reported 
range  of  f(u)  (Fig.  14)  is  insufficient  to 
verify  where  this  happens  for  the  lower 
E/N  values.  Other  electron-induced  rates, 
such  as  vibrational  excitation  and 
dissociation,  must  be  also  strongly  affected 
by  the  distortion  in  f(u). 

The  traditional  way  to  account  for 
vibrational  (more  generally,  inelastic) 
losses  in  molecular  gases  is  to  multiply  the 
elastic  losses  times  a  factor  5,  which  is 
energy-dependent.  Data  on  these  factors 
were  collected  by  Massey  and  Craggs 
and  are  reported  in  Ref.  (32),  for  example. 
Eq.  (5.74)  of  Ref.  32  gives  the  electron 
temperature  elevation  as 


T,m,  AfY 

2AST.m,Q„(T.){F) 


(15) 


We  can  use  the  Nighan  results  in 
Fig.  (14)  and  Eq.  (15)  (w'here  P=N  k  Tg)  to 
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extract  the  effective  6  value  implied.  We 
do  this  by  equating  Te  to  the  mean  energy 

u,  in  Fig.  (14),  and  using  known  data  for  the 
elastic  cross-section  Qeg  (Tg)  (Ref.  45).  The 
value  of  Tg  needs  to  be  assumed,  and  0.1  eV 
has  been  used  here.  The  resulting  6  is 
shown  superimposed  on  the  Massey-Craggs 
data  in  Fig.  16,  and  a  very  large 
discrepancy  is  readily  appreciated.  With 
Nighan's  results,  if  E/N-IO*^^  Vcm^,  as  in 
the  anodic  bridge  region,  then  Te-0.6  eV, 
and  a  very  non-maxwellian  distribution 
results,  while,  using,  as  in  Refs.  (20-31),  the 
Massey-Craggs  data,  we  found  Te~l-2  eV, 
which,  back  to  Fig.  14,  would  indicate 
little  departure  from  Maxwellian. 


Fig.  16:  Energy  loss  of  electrons  with 
maxwellian  velocity  distribution  with 
various  gases.  (Data  compiled  by  H. 
Massey  and  J.D.  Craggs  (44)  1959). 
Superimposed  are  calculations  on  8  based 
on  results  in  Fig.  14  for  N2. 


In  view  of  the  uncertainties,  and  of 
the  strong  potential  implications,  there  is 
clearly  a  need  for  more  detailed  research  in 


this  particular  area.  As  noted,  Babu  and 
Subramaniam  have  initiated  work  in 
that  direction. 

iA  Thermo-Mechanical  Effects 

The  materials  comprising  the 
anode  block  structure  of  an  arcjel  are  subject 
to  extreme  heat  fluxes  (1-lOxiO^  w/mMn 
the  attachment  area),  and  reach 
temperatures  that  can  be  over  1/2  of  their 
melting  point.  Thus,  despite  a  low  stress 
level,  thermal  creep  effects  can  occur, 
which  distort  the  critical  constrictor 
region,  as  has  been  observed,  for  example  by 
Butler  and  lead  to  performance  drift 
and  lack  of  repeatability,  and  may  also 
affect  thruster  life.  Related  effects  can 
occur  in  cathodes  at  high  current  From 
the  modeling  point  of  view,  these  areas 
remain  unexplored,  and,  given  their 
practical  importance,  efforts  should  be 
made  to  extract  as  much  related 
information  as  possible  from  heat  flow 
models  of  the  solid  part  of  thrusters. 

4.5  Species  Segregation  and  Mean  Free 
Path  Effects 

In  an  arejet  using  N2H4  or  NH3, 
these  heavy  molecules  will  dissociate  as 
they  enter  the  arc-heated  parts  of  the 
flow.  The  hydrogenic  fragments  resulting 
from  this  dissociation  will  diffuse  much 
faster  than  the  nitrogen  molecules  or  atoms, 
and  the  result  can  be  a  hydrogen  enrichment 
of  the  cooler,  outer  parts  of  the  flow. 

A  similar  segregation  of  hydrogen 
towards  the  periphery  is  also  to  be 
expected  in  the  plume  after  collisionality 
has  become  very  weak,  simply  because  of 
the  higher  thermal  velocity  of  H  and  H2 
compared  to  N  and  N2  (or  NH)  at  the  same 
temperature. 

These  two  effects  are  observed 
together  when  looking  spectroscopically  at 
the  arejet  plume,  yet  their  implications  are 
quite  distinct  (segregation  upstream  can 
improve  specific  impulse,  segregation  in 
the  plume  cannot).  Refined  diffusional 
and/or  Monte  Carlo  computations  will  be 
required  to  mimic  these  effects. 
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5.  Some  Numcrica]  Tssup.; 

With  so  many  difficult  physical 
problems  pending,  arcjet  modelers  have  not 
so  far  paid  much  attention  to  questions  of 
numerical  efficiency.  The  standard 
approach  for  2-D  simulations  has  been  to 
combine  an  explicit  time-marching 
algorithm,  such  as  McCormack’s,  for  the 
bulk  flow,  plus  over-relaxation  or  other 
elliptic  solvers  for  the  field  equations. 
Miller  found  it  necessary  to  integrate 
the  electron  continuity  and  energy  equations 
and  the  potential  equation  using  a  separate 
time  step,  which  needs  to  be  up  to  100  times 
shorter  than  that  used  for  the  heavy- 
particle  equations.  The  time  steps  used 
need  to  satisfy  the  Courant-Fredrich-Levy 
condition  (the  fastest  signal  should  travel 
less  than  one  mesh  size  per  time  step).  In  a 
very  inhomogeneous  situation,  such  as  an 
arcjet,  the  fastest  velocities  and  speeds  of 
sound,  which  occur  in  the  hotter  parts  of 
the  nozzle  flow,  dictate  this  choice  of  time 
step,  with  the  result  that  the  slow  inlet 
region  is  enormously  over-calculated. 
Combining  these  inefficiencies  with  the 
abundance  of  exponentials  and  other 
computationally  costly  calculations 
required  per  step,  the  most  advanced  codes 
are  at  present  pushing  against  the  frontiers 
of  what  is  practical  in  terms  of  computation 
costs. 

It  seems  clear,  therefore,  that  the 
time  has  come  to  begin  refining  the 
computational  procedures,  if  we  want  to 
continue  enriching  the  codes  with 
additional  effects,  species  and,  especially, 
dimensions.  This  is  definitely  beyond  this 
writer’s  expertise,  but  a  few  indications  as 
to  improvements  to  be  sought  will  be 
ventured  anyway: 

(a)  Simplify  and  pre-tabulate  as 
many  of  the  transport  and  rate  coefficients 
as  possible. 

(b)  Try  to  incorporate 
inhorrtogeneous  time  step  procedures,  if  only 
the  steady  state  is  of  interest. 


(c)  Use  variations  of  the 
unstructured  grid  algorithms,  now  being 
developed  for  CFD  applications. 

(d)  Investigate  the  use  of  quasi¬ 

equilibrium  for  the  fastest  evolving  effects, 
such  as  electron  temperature 

balance. 

6.  Summary 

Numerical  modeling  of  arcjet 
thrusters  is  reaching  an  acceptable  level  of 
performance  prediction  capacity, 
predicated,  however,  n  judicious  admixing 
of  empirical  information  on  a  number  of 
poorly  understood  physical  issues. 
Prominent  among  these  are  the  physics  of 
anodic  attachment  (including  stability  and 
potential  filamentation),  and  the 
interactions  of  free  electrons  with  the 
molecular  gas.  Further  progress  will 
necessarily  have  to  include  clarification  of 
these  areas,  plus  gradual  expansion  of  the 
ability  to  handle  the  very  complex 
numerical  task  while  staying  financially 
solvent.  The  time  appears  at  hand  when 
theoretical  modeling  can  emerge  from  the 
shadow  of  empiricism  and  become  a  co¬ 
equal  driver  with  selective  testing  in  the 
quest  for  improved  design  of  future 
thrusters. 
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Quantity 

Computational 

(Experimental) 

m(g/scc) 

0.1 

0.1 

0.3 

0.3 

1  (A) 

50 

150 

50 

150 

Po(N/m2) 

9.02  xlCH 

1.23x105 

2.03x1^ 

2.44X105 

(9.00  xl0<) 

(1.30x105) 

(2.35  xl05) 

(2.74X105) 

V  (V) 

165 

155 

225 

199 

(165) 

(140) 

(233) 

(198) 

Thrust  F  (N) 

OSS 

0.81 

132 

1.63 

(0.60) 

(0.83) 

(1.17) 

(1.62) 

TJ 

0.204 

0.141 

0358 

0.148 

(0.218) 

(0.164) 

(0.196) 

(0.147) 

Qloss  (KW) 

1.40 

4.20 

1.40 

(1.30) 

(4.00) 

(1.00) 

Table  1.  Comparison  of  1-D  model  to  data  H^for  a  20  KW  water-cooled  H2  arcjet.  A 
uniform  anode  drop  of  AV,  =  28^.  has  been  added  to  the  calculated  column  voltage. 
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Abstract 

A  detailed  numerical  model  has  been  developed  to 
study  the  gasdynamic  flow  in  an  electrothermal  m- 
cjet  thruster.  This  twc^temperature,  Navier-Stokes 
model  consistently  incorporates  viscosity,  heat  con¬ 
duction,  ohmic  dissipation,  coUisional  energy  trans¬ 
fer  between  electrons  and  heavy  species,  ambipo- 
lar  diffusion,  nonequiUbrium  dissociation  and  mira¬ 
tion,  and  radiation.  The  fluid  equations  are  solved 
by  MacCormack’s  method,  while  an  iterar®  P*®" 
cedure  is  used  to  relax  an  electric  potential  equ^ 
tion,  from  which  the  current  distribution  in  the 
thruster  is  obtained.  Converged  solutions  a^  com¬ 
pared  with  experimental  results  &om  the  G«man 
TTl  radiativcly-cooled  arcjet  thruster  with  hydro- 
gen  propeUant.  Results  are  presented  for  a  bas^ 
line  case  which  reveal  the  twc^dimensiond,  twc^fluid 
nature  of  the  interior  flow,  especially  m  terms  of 
the  distribution  and  anode  attachment  of  the  ele^ 
trie  current  and  the  growth  and  development  of  the 
arc  region.  Calculated  discharge  voltage  is  within  a 
few  percent  of  experimental  measurements,  and  pre¬ 
dicted  specific  impulse  is  within  5*10%  ag^emen 
over  a  range  of  operating  parameters.  The  effects  of 
a  coupled  anode  heat  balance  model  on  the  predicted 
anode  wall  temperature,  inlet  gas  temperature,  and 
overall  thruster  performance  are  also  discussed. 
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Vibrational  excitation  energy  [J] 
Planck’s  constant  [J<],  enthalpy  [J/kfl] 
Total  current  [A] 

Electric  current  density  [Airr?] 
Boltzman’s  constant  [J I^K\ 


Equilibrium  constant 
Particle  mass  [kp] 

Number  density  [1/ru®] 

Net  production  rate  [l/m*/s] 
Avogadro’s  number  [l/mole] 

Scalar  pressure  [Po] 

Heat  flux  vector  [PV/m*] 

Real  gas  constant  (J/kp/*if] 
Universal  gas  constant  [J ImcAcj  K] 
Energy  loss  due  to  radiation  [kV/m  ] 
Temperature  [®X] 

Mean  flow  velocity  [m/s] 

Bohm  velocity  [m/s] 

Slip  velocity  [m/s] 

Mole  fraction 
loniz&tion  fraction 


Coefficient  of  viscosity  [kp/ms] 

Collision  &equency  [1/s] 

Electric  potential  [V]  ^ 

Viscous  dissipation  function  [W/m®] 
Electron  mobility  [m*/0/C] 

Mass  density  [kp/m®] 

Electrical  conductivity[mho/m] 
Average  effective  coUision  integral  [m  ] 


1  Introduction 

Low  power  arcjet  thrusters  have  recently  been  flight 
qualified  through  ground  testing  and  will  soon  be 
tested  in  space  for  stationkeeping  of  geosynchronous 
communication  sateUites.  Most  of  the  impetus  for 
design  strategies,  however,  has  come  from  empirical 
studies  and  experimentation,  and  a  need  remains  to 
better  understand  the  underlying  phyrics,  detailed 
energy  balances,  and  transport  mechanisms  of  these 
devices.  Unfortunately,  the  complexity  of  the  mod¬ 
els  and  equations  needed  to  accurately  represent  the 
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flow  in  an  aicjet  thruster  effectively  limits  the  use 
of  analytic  techniques  to  simplified  cases,  through 
which  one  may  obtain  useful  physical  insights  but 
inadequate  predictions  of  thruster  performance.  Ex¬ 
perimental  techniques  provide  much  useful  empirical 
data,  but  many  quantities  of  interest  are  not  acces¬ 
sible  b  the  important  regions  of  the  thrusters.  For 
these  reasons  numerical  methods  of  solvbg  the  gov- 
ernbg  equations  have  become  an  important  tool  for 
conduetbg  arejet  research. 

Previous  numerical  modelbg  of  arejet  thrusters 
has  focused  on  the  development  of  1-D,  2-D,  and 
axisymmetric  models  with  relatively  simple  physics 
and  geometries.  The  level  of  detail  has  ranged 
&om  1-D  models[l]  to  coupled  quasi-analytic  mod¬ 
els  of  the  bner  (arc)  and  outer  (cold  gas)  flows[2, 
3,  4],  to  simplified  axisymmetric  space-marchbg 
techniques [5],  and  finally  to  2-D  and  axisymmet¬ 
ric  viscous  codes  which  begb  to  bcorporate  most 
of  the  detailed  physical  processes[6,  7,  8],  The  lat¬ 
est  research  has  obtabed  results  which  variously  b- 
cludc  ohmic  heatbg,  electron  heat  conduction,  ther¬ 
mal  and  ionbational  noncquilibrium,  and  empirical 
models  of  radiation  losses.  There  are  still  a  num¬ 
ber  of  issues,  such  as  viscous  and  diffusive  effects, 
arc  formation  and  attachment,  the  heat  balance  b 
the  anode,  and  ultimately  the  accurate  prediction  of 
voltage  and  efficiency,  which  need  to  be  addressed. 
This  paper  describes  a  generalized,  more  physically 
accurate  model  of  the  gasdynamic  flow  through  an 
arejet  thruster,  which  bcludes  the  aforementioned 
effects  and  compares  favorably  to  experimental  re¬ 
sults  obtabed  with  medium  power  hydrogen  arejets. 


2  Model 

2.1  Basic  Assumptions 

This  model  is  based  on  an  axisymmetric  formulation, 
so  that  variations  b  flow  quantities  b  the  azimuthal 
direction  are  neglected.  A  component  of  the  flow  ve¬ 
locity  b  the  ^-dbection,  however,  is  bcorporated  to 
account  for  the  “swbF  bjection  of  most  experimen¬ 
tal  arejets.  The  model  has  been  developed  b  a  gen¬ 
eral  enough  sense  so  that  any  monatomic  or  diatomic 
propellant  may  be  simulated.  For  the  purpose  of 
this  research  hydrogen  was  selected  as  the  propel¬ 
lant  of  choice,  due  to  its  low  molecular  weight  (and 
therefore  high  performance)  and  its  simple  molecu¬ 
lar  structure,  which  allows  for  analytic  evaluation  of 
the  necessary  transport  coefficients.  Nonequbbrium 
dissociation  and  ionization  are  modeled,  and  four 
species  of  particles  are  tracked:  diatomic  molecules, 
monatomic  neutrals  and  ions,  and  electrons.  Disso¬ 
ciation  is  modeled  by  heavy  species  collisions  and  by 
electron  impact,  and  the  ionization  process  is  based 


on  electron  impact  ionbation  and  three-body  recom- 
bbation,  with  only  H'^  ions  considered. 

The  foUowbg  assumptions  are  made  regarding  the 
state  of  the  flow  b  the  thruster  and  the  physical 
processes  bvolved.  The  plasma  produced  by  ioniz- 
bg  electron  collisions  is  assumed  to  be  macroscop- 
ically  neutral,  so  that  n«  =  Strong  couplbg 
is  assumed  between  the  ions  and  neutrals,  desig¬ 
nated  together  as  the  heavy  species.  This  implies 
that  tZi  =  tin  2?  tT  (except  for  ambipolar  diffusion), 
and  Ti  ^  Tn  ^  Effects  which  arc  consistently 
bcorporated  bclude  ambipolar  diffusion,  heat  con¬ 
duction,  viscous  shear  and  dissipation,  ohmic  heat- 
bg,  coUisiorval  energy  transfer  between  electrons  and 
heavy  species,  and  energy  lost  through  radiation. 
The  self-bduced  magnetic  field  of  the  ionbed  gas 
b  neglected  due  to  the  low  current  density  b  the 
thruster,  and  the  bdividual  species  arc  assumed  to 
obey  the  ideal  gas  law.  Given  the  aforementioned 
assumptions,  the  model  can  be  summarbed  by  a  set 
of  nbe  partial  differential  equations  which  must  be 
solved  locally  b  order  to  generate  a  viable  simula¬ 
tion  of  the  flow  b  an  arejet  thruster. 


2.2  Governing  Equations 

The  set  of  equations  which  govern  the  flow  b  the 
model  arejet  thruster  of  thb  research  is  essentially 
a  group  of  modified  Navier- Stokes  equations.  These 
bdude  an  equation  of  state  and  equations  for  the 
ion,  neutral  atom,  and  global  density;  the  axial,  ra¬ 
dial,  and  azimuthal  global  momentum;  the  electron 
and  heavy  species  energy;  and  the  electric  potential. 
In  order  to  mabtab  numerical  robustness,  these 
equations  are  written  b  as  conservative  a  manner 
as  possible,  particularly  with  respect  to  the  ~  terms 
which  appear  due  to  the  axbymmetry  of  the  prob¬ 
lem. 

The  electric  potential  equation  is  derived  by  com- 
bbbg  Ohm’s  law,  J  =  with  the  equation 

V  •  j  =  0,  where  ^  =  "■'»  * - b  the  electron  mo- 

biiity.  Assumbg  a  potential  of  the  form  E  =  — 
the  resultbg  equation  b  given  by 
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~  ^  I  fir*  r  dr  dz^  )  dr  dr  dz  dz 

(1) 

The  current  density  may  then  be  extracted  by  solv¬ 
ing 


Jr  =  V- 


dpt 

dr 


and 


Jz 


(2) 


2 


The  global  density  equation  b  obtained  by  sum¬ 
ming  the  individual  species  continuity  equations: 

dpr  dpUrT  8pu,r  _ 

&t  ^  8r  dz 

For  the  ion  (or  electron)  density,  the  governing  equa¬ 
tion  b  derived  from  the  statement  of  ion  mass  con¬ 
servation,  modified  to  account  for  the  ambipolar  flux 
of  charged  particles  in  terms  of  ion  density  gradients: 

SptT  d(p,Ur  +  dtt)r  ^  d{p,u,+d,,)r  _ 

^  (4) 

where  d,  b  the  ambipolar  flux,  given  by 

d.  -  -DcVp,-,  -  Y  4^  +  nn)’ 

(5) 

Here  the  source  term  n,  represents  the  net  rate  of 
production  of  ions  per  unit  volume  through  inelas¬ 
tic  coUbional  processes.  The  statement  of  atomic 
hydrogen  mass  conservation  b  obtained  in  a  simi¬ 
lar  manner,  and  the  effect  of  ambipolar  diffusion  b 
consbtently  included  by  assuming  that  both  neutral 
species  travel  at  the  same  velocity: 

flpjg-r  SpBUfT  _  ^  (  PB 

8r  8r  \{pb3+Pb)  "  / 

8pBV,r  f  PB  ^B 

8z  8z  \ {pB,  +  Pe)  ) 

—  rriH  (fisd"  ^  ^  ~  fi«)  (®) 

In  thb  equation  the  source  term  represents  the  net 
rate  of  production  of  neutral  atoms,  given  by  the 
difference  between  the  net  rate  of  production  due  to 
dbsociation  and  the  net  rate  of  ionbation  of  those 
neutral  atoms  produced. 

The  three  global  momentum  equations  are 

6pUrT  e(pu*  -l-p-  T„)r  _  g(p«TUx  -Trx)^’ 
__  +  _  + 


+  -  P  =  0, 

'pUrUt  —  Trt)T  8[f»t, 

8r 

-(-  pUrU$  =  0,  and 


8fm$r  d{pUrU0 -- Tfe)r  ^  djpu^u^  —  T(p^)r 

ez 


dpu^r  8(pUrUt  —  Vti)t  8(pul  +  p-  r»f)r  _  ^ 
dr  8z 

where  the  corresponding  are  the  viscous  stress 
tensor  components  in  axisymmetric  coordinates. 
The  heavy  species  energy  equation  is 

8pt„r  8{pesUr  +qar)r  8[pe,u,  +q,t)r 

-JT  ^  Tr  ^  8z 


8  {pBi^Bj^Btr  •+  PB^B^Bt  PE*^B*^B*r)^ 
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where  the  internal  energy  and  enthalpy  are  defined 
by 

pe,  =  ^{pB  +  PB*)  ^B^f  +  :^PE3RB3Tg 

(11) 
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phg  =  p€g  -f  Pff^Rff^Tg  +  {PB  ^  ^B^g »  (l2) 

Additional  flux  terms  have  been  included  in  the  ra¬ 
dial  direction  to  account  for  the  radial  transport  of 
energy  by  ambipolar  diffusion.  The  heat  flux  vectors 
in  Eqn.  10  can  be  expressed  as 

=  =  (IS) 

while  the  viscous  dissipation  function  is  given  by 
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Collisional  energy  transfer  between  electrons  and 
heavy  species  is  represented  by  the  following  expres¬ 
sion: 

El  =  -h  UtH  +  ^ 

ms 

The  coefficient  6^  in  Eqn.  15  is  necessary  to  correct 
for  the  fact  that  electron-F:  collisions  are  inelastic 
in  nature. 

For  electrons,  the  governing  energy  equation  is: 

8p,E,r  8{p,u,rBt  +  qtr)T  8(p,VttSt  +  q,z)r 

dt  8r  8z 


—  fi —  El  -  Ei  <  o-v  >  n.njj,  -  Eih,  -  .rJ  r. 

(16) 
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Here  ^  represents  heating  due  to  ohmic  dissipation, 
R  is  the  contbuum  radiative  loss,  and  the  total  en¬ 
ergy  and  enthalpy  are  given  by 

i:<  =  Ir.Z  +  \uI  (17) 

and 

Re=:£,^R,Te.  (18) 

The  required  electron  velocities  are  extracted  from 
the  local  current  density. 

The  final  equation  required  for  closure  of  the  set 
is  the  equation  of  state: 

p  =  -f  n<iTe . 

J  3 

2.3  Dissociation  and  Ionization 

The  nonequilibrium  dissociation  rate  ns  is  derived 
by  following  the  procedure  and  nomenclature  of 
Biasca[9]  for  collisions  of  heavy  species.  Accordingly, 
this  rate  is  given  by 

tie  =  ANT^  eip  j  +  mE,<rE^) 


X 


RT, 

K,(T,) 


(20) 


where  N  is  Avogadro’s  number,  Kj,  is  the  equilib¬ 
rium  constant  in  terms  of  partial  pressures,  the  CjS 
are  the  species  molar  concentrations,  and  the  ap¬ 
propriate  constants  for  hydrogen  are  listed  in  Ta¬ 
ble  1[10].  Dissociation  by  c  —  fi’2  collisions  is  repre¬ 
sented  by  the  second  term  on  the  right-hand  side  of 
Eqn.  6,  where  the  reaction  rate  coefficient  <  (rv  > 
as  a  function  of  Tt  is  taken  from  Janev  et  ed.[Il]. 


Table  1:  Constants  for  the  Hydrogen  Dissociation 
Rate  Equation _ 


Constant 

Value 

A  {m^/mole  —  s) 

5.5  X  10^^  ' 

B  (J/mole) 

435,600 

n 

-1 

^E 

5 

rhE, 

2 

For  ionfration,  the  finite  production  rate  is 
given  by  the  generalked  model  of  ionization 
and  three-body  recombination[12]  as  modified  by 
Sheppard[l3]: 


n,  =  Rut  (Sub  -n*) 


(21) 


S  = 


f  2xTn 


V 


R  =  6.985  X 


13^  -  4-0833)^ 
0.8179 


(23) 


2.4  Transport  Properties 

Because  of  the  multi-component  nature  of  the  gas 
in  an  electrothermal  arejet,  even  with  hydrogen  as 
a  propellant,  the  equations  for  the  transport  prop¬ 
erties  become  quite  complex.  Since  data  are  more 
readily  available  for  hydrogen  in  the  form  of  colli¬ 
sion  integrals,  a  formulation  of  the  transport  coeffi¬ 
cients  based  on  these  integrals  rather  than  on  colli¬ 
sion  cross-sections  is  implemented.  As  an  example 
of  the  form  of  the  transport  coefficients  used  in  this 
model,  the  gas  coefficient  of  viscosity  is  given  by 
Eqn.  24.  It  is  calculated  based  on  a  mean  free  path 
mixture  formula,  which  is  a  function  of  the  collision 
integrals  and  the  species  number  densities  and  pure 
viscosities[12]: 
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where  the  pure  viscosities  are  given  to  lowest  order 
by 

Hi  =  2.6693  X  10~“  .  (25) 

(y  ■  ’ 

**33 

Electron  and  heavy  species  thermal  conductivity  co¬ 
efficients  are  calculated  based  on  similar  mean  free 
path  arguments  and  mixture  rules,  and  the  electrical 
conductivity  is  taken  from  a  first  order  approxima¬ 
tion  by  Grier [14].  Collision  integrals  required  in  the 
calculation  of  transport  coefficients  are  interpolated 
as  a  function  of  temperature  from  data  by  Grier  [14], 
Belov[15],  and  Vanderslice  ct  al.[l6].  The  accuracy 
of  the  above  approximate  formulas  for  the  transport 
coefficients  was  verified  by  comparison  to  previous, 
more  detailed  work  for  hydrogen  in  thermal  and  ion- 
izational  equilibrium. 


3  Numerical  Method 

3.1  Integration  Scheme 


(22) 
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The  govetniiig  fluid  equations  are  numerically  in¬ 
tegrated  using  MacCormack’s  explicit  node-based 
method,  which  is  2nd  order  accurate  in  time  and 


space.  An  empirical  stability  formula  developed  by 
MacCormack  and  Baldwin[17]  is  utiliied  to  calculate 
the  integration  time  step,  which  is  then  multiplied 
by  a  fractional  coefficient  to  account  for  the  effect  of 
source  terms  and  nonlinearities.  Although  MacCot- 
mack’s  method  contains  some  inherent  dissipation, 
additional  numerical  smoothing  is  requited  in  order 
to  damp  unwanted  numerical  oscillations  in  the  fluid 
equations.  Consequently,  2nd  and  4th  order  smooth¬ 
ing  terms  ate  applied  to  each  equation  as  necessary. 
The  electric  potential  equation,  being  predominantly 
elliptic  in  nature,  is  solved  by  iteration  using  a  suc¬ 
cessive  overrelaiation  (SOR)  technique. 

Since  the  physical  grid  is  nonuniform  in  order  to 
closely  represent  the  geometry  of  the  actual  arcjet 
being  modeled,  the  governing  equations  are  trans¬ 
formed  into  natural  coordinates  and  then  solved  on 
a  uniform  Cartesian  computational  mesh.  The  phys¬ 
ical  grid  contains  140  axial  by  30  radial  grid  points. 

3.2  Boundary  Conditions 

The  conditions  at  the  inlet  of  the  computational  do¬ 
main  are  postulated  to  be  those  of  a  flow  which  has 
just  been  injected  into  the  thruster  plenum  by  a  large 
number  of  evenly  spaced  jets.  The  flow  is  therefore 
assumed  to  be  subsonic  and  parallel  to  the  thruster 
walls.  A  fraction  of  the  total  inlet  velocity,  typically 
30-50%,  is  specified  as  being  in  the  aiimuthal  di¬ 
rection  in  order  to  simulate  an  injected  swirl.  The 
mass  flow  rate  and  total  enthalpy  are  specified  based 
on  the  particular  run  parameters,  and  the  ioniiation 
fraction  is  set  to  a  small  value,  typicaUy  1  x  lO"*. 
The  density  is  obtained  from  a  downwind  finite  dif¬ 
ference  approximation  of  the  overall  continuity  equa¬ 
tion,  and  the  inlet  electron  temperature  is  set  equal 
to  that  of  the  next  inside  point.  No  current  is  al¬ 
lowed  to  pass  upstream  of  the  inlet. 

The  boundary  conditions  at  the  outlet  of  the 
thruster  depend  on  whether  the  exit  flow  is  subsonic 
or  supersonic.  In  both  cases  the  electron  tempera¬ 
ture  is  set  equal  to  that  of  the  next  i^ide  point,  and 
no  current  may  pass  beyond  the  exit  plane.  If  the 
flow  is  supersonic  at  a  point  on  the  exit  plane,  then 
the  remaining  quantities  are  extrapolated  from  their 
values  at  the  preceding  two  grid  points  of  the  mesh. 
If  the  flow  is  subsonic,  then  the  exit  pressure  is  set 
equal  to  a  small  value  representing  near-vacuum  con¬ 
ditions  and  the  density  and  axial  velocity  are  given 
by  Riemann  invariants.  The  remaining  quantities 
at  a  subsonic  outlet  are  then  calculated  as  in  the 
supersonic  case. 

For  those  boundary  points  lying  beyond  the  tip  of 
the  cathode  on  the  line  of  symmetry  (r  =  0),  the 
radial  and  azimuthal  flow  velocities  are  set  equal  to 
zero  and  a  zero  radial  gradient  is  imposed  on  the 


remaining  quantities. 

Flow  boundary  conditions  at  the  thruster  walls 
include  viscous  no-slip  conditions  for  the  fluid  ve¬ 
locities  and  an  imposed  zero  gradient  on  the  elec¬ 
tron  temperature.  The  heavy  species  temperature  is 
held  constant  at  1000®K  upstream  of  the  constric¬ 
tor,  increasing  linearly  to  1100®K  at  the  constric¬ 
tor  exit.  This  profile  was  chosen  based  on  experi¬ 
mental  and  numerical  calculations  of  the  anode  wall 
temperature  distribution  for  a  reasonable  operating 
range  of  the  German  TTl  radiatively-cooled  arcjet 
thruster [18].  On  the  cathode  the  wall  temperature  is 
allowed  to  increase  to  a  maximum  of  2000‘’K  at  the 
tip.  For  the  boundary  condition  on  electron  density 
at  each  electrode,  a  balance  is  postulated  between 
the  flux  of  ions  arriving  at  the  sheath  edge  by  am- 
bipolar  diffusion  and  the  flux  of  ions  arriving  at  the 
wall  by  virtue  of  their  thermal  energy  at  the  Bohm 
velocity  (tin): 

=  O.eintVg,  (26) 

dy 


where 


This  boundary  condition  neglects  the  voltage  drops 
present  in  the  non-neutral  plasma  sheath.  Also 
omitted  is  the  dependence  of  the  thermionic  emis¬ 
sion  of  electrons  at  the  cathode  on  the  cathode  wall 
temperature.  Application  of  the  perpendicular  over¬ 
all  momentum  equation  at  the  walls  provides  the 
approsmate  condition 

^|«aZI=0,  (28) 

an 

where  inertial  and  viscous  terms  have  been  ne¬ 
glected. 

The  boundary  condition  on  the  wall  electric  po¬ 
tential  is  that  there  is  no  current  perpendicular  to  an 
insulating  section,  and  that  the  potential  on  the  an¬ 
ode  is  set  equal  to  a  fixed  but  arbitrary  voltage.  For 
numerical  reasons,  anode  current  attachment  is  re¬ 
stricted  to  that  portion  of  the  outer  wall  downstream 
of  the  constrictor  exit.  On  the  cathode  tip,  a  uniform 
nrifil  current  density  is  prescribed  which  sums  to  the 
specified  total  current;  the  potential  at  the  cathode 
is  then  chosen  so  as  to  maintain  this  current  level. 
A  cathode  voltage  drop  equal  to  the  ionization  po¬ 
tential  plus  one  half  of  the  dissociation  potential  of 
the  gas  is  added  to  the  calculated  voltage  in  order 
to  account  for  the  model’s  neglect  of  the  cathode 
sheath  region.  In  addition,  an  anode  voltage  drop 
is  subtracted  from  the  calculated  voltage  in  order  to 
account  for  the  anode  sheath  region.  A  negative  po¬ 
tential  gradient  is  required  in  this  sheath  in  order  to 
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turn  back  excess  electrons  since  the  extracted  cur¬ 
rent  is  much  less  than  the  random  thermal  flux  of 
electrons  to  the  anode  wall  (j»nedt  <  ^en,c,).  This 
voltage  drop  is  given  by 


4  Results 

4.1  Baseline  Case 

Results  have  been  achieved  foi  comparison  to  the 
German  TTl  radiation-cooled  arcjet  thruster[l8]  at 
a  number  of  operating  points.  For  the  baseline  case 
of  J  =  100>1  and  m  =  O.ly/s,  some  results  have 
been  presented  in  a  previous  paper [19].  Table  2  com¬ 
pares  bulk  results  of  the  simulation  to  experimen¬ 
tal  measurements,  while  Figures  1  through  9  show 
line  and  contour  plots  of  representative  quantities. 
There  is  exceUent  agreement  in  the  voltage  between 
predicted  and  experimental  results,  and  good  agree¬ 
ment  (within  7%)  in  the  thrust  and  specific  impulse. 
The  accuracy  in  voltage  prediction  shows  that  the 
model  is  accurately  modeling  arc  growth,  electrical 
conductivity,  and  current  attachment  given  the  as¬ 
sumptions  we  have  made  regarding  electrode  sheath 
voltage  drops.  The  overprediction  of  thrust  may  be  a 
result  of  inaccurately  specified  boundary  conditions 
for  the  electrode  temperatures  and  the  inlet  gas  tem¬ 
perature  (see  Sec.  5),  or  it  may  reflect  remaining 
deficiencies  in  modeling  the  anodic  arc  attachment 
region. 


Table  2:  Comparison  of  Predicted  and  Experimental 
Results  for  Baseline  Case  _ _ 


Predicted 

Experiment 

Voltage  (V) 

115 

112 

Power  (kW) 

11.5 

11.2 

Thrust  (N) 

1.01 

0.94 

Specific  Impulse  (s) 

1030 

960 

Efficiency 

0.442 

0.395 

Figure  1  shows  an  axial  line  plot  of  the  electric 
potential  from  the  cathode  tip  to  the  anode  attach¬ 
ment  rone.  The  near-cathode  voltage  drop  AV^  is 
composed  of  a  15.8V  drop  assigned  to  the  cathode 
sheath  (V-  ^  \Vi)  and  an  8V  drop  calculated  in  the 
first  few  grid  points  downstream  of  the  tip.  The 
near-anode  voltage  drop  AK  %  15V  is  composed  of 
a  22V  drop  captured  by  the  simulation  and  a  -7V 
drop  associated  with  the  electron-repelling  sheath. 


This  total  anode  voltage  drop  is  associated  with  the 
net  deposition  of  energy  into  the  anode  block  by 
heavy  species  heat  conduction  and  by  the  impinge¬ 
ment  of  current-carrying  electrons.  Assuming  that 
the  energy  transferred  per  unit  area  is  of  the  form 

-f- , 

(30) 

using  the  results  of  the  baseline  flow  simulation 
yields  an  equivalent  voltage  of  14.5V  for  this  de¬ 
posited  power,  which  agrees  weD  with  the  AKi  seen 
in  the  potential  profile. 

Current  streamlines  are  plotted  in  Figure  2.  In^ 
this  case,  the  bulk  of  the  current  attaches  within  the 
first  quarter  of  the  nozzle,  with  a  peak  just  down¬ 
stream  of  the  constrictor  exit.  The  flow  becomes 
fully  ionized  along  the  centerline  immediately  down¬ 
stream  of  the  cathode  tip  and  remains  so  through 
the  first  part  of  the  nozzle  expansion,  beyond  which 
there  is  some  recombination  (Figure  3).  The  bound¬ 
ary  of  the  partially  ionized  region  grows  to  approx¬ 
imately  50%  of  the  channel  by  the  constrictor  exit, 
and  this  region  is  essentially  entrained  in  the  flow 
throughout  the  nozzle.  The  primary  heating  mecha¬ 
nism  is  ohmic  dissipation,  which  peaks  locally  along 
the  constrictor  centerline  and  just  beyond  the  con¬ 
strictor  exit  near  the  anode.  This  is  evidenced  by 
the  local  maxima  in  electron  temperature  in  these 
regions  (Figure  4). 

Within  the  highly  ionized  region  of  the  arc  in  the 
constrictor,  collisional  energy  transfer  between  elec¬ 
trons  and  heavy  species  raises  the  gas  temperature 
to  20, 000  -  30,  OOO^iT,  or  nearly  the  same  temper¬ 
ature  as  the  electrons.  Outside  of  the  arc  the  flow 
remains  cold,  at  a  temperature  approximately  equal 
to  the  anode  wall  temperature.  This  outer  flow  re¬ 
mains  essentially  uncoupled  from  the  hot  core  flow, 
as  pictured  in  Figure  5.  The  flow  velocity  is  also  un¬ 
coupled,  although  viscous  forces  eventually  wash  out 
the  separation  in  the  nozzle  expansion  (Figure  6). 
Since  the  pressure  is  nearly  uniform  in  the  radial 
direction,  rapid  acceleration  of  the  core  flow  occurs 
throughout  the  low  pressure  region  of  the  arc  in  the 
constrictor.  Once  the  bulk  of  the  pressure  work  has 
been  utilized  in  the  expansion  process,  however,  vis¬ 
cous  forces  arising  &om  steep  velocity  gradients  in 
the  central  core  decelerate  the  flow  significantly  in 
the  nozzle  expansion.  Both  the  inner  and  outer  flows 
accelerate  smoothly  through  sonic  velocity  at  and 
just  beyond  the  constrictor  exit,  and  a  peak  Mach 
number  of  2.86  is  reached. 

Simulating  the  current  attachment  at  the  anode 
realistically  and  self-consistently  has  been  a  major 
difficulty  in  previous  arcjet  simulations.  The  effec¬ 
tiveness  of  this  model  in  simulating  this  region  is 
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dne  to  the  incorporation  of  separate  energy  equa¬ 
tions  for  the  heavy  species  and  electrons  and  to  the 
nse  of  nonequilibrium  dissociation  and  ionization  fi¬ 
nite  rate  equations.  Ohmic  dissipation  is  found  to 
be  an  important  source  of  energy  both  inside  and 
outside  the  arc  (Figure  7).  This  leads  to  electron 
temperatures  as  high  as  20,000'’if  in  the  anode  at¬ 
tachment  zone  of  the  outer  flow,  much  higher  than 
the  1000-2000®^r  temperatures  which  would  be  cal¬ 
culated  by  a  model  with  only  one  energy  equation. 
Figure  8  illustrates  this  result  by  comparing  radial 
profiles  of  the  electron  and  heavy  species  tempera¬ 
tures  at  an  axial  location  0.25mm  downstream  of  the 
constrictor  exit.  This  elevated  electron  temperature 
then  produces  enough  electron  impact  dissociation 
and  ionization  to  create  the  necessary  charge  carri¬ 
ers  for  electrical  conduction  between  the  outer  arc 
boundary  and  the  anode.  Radial  ambipolar  diffu¬ 
sion  also  plays  a  role  in  moving  ions  and  electrons 
outward  &om  the  arc  into  the  surrounding  cooler  gas 
flow.  This  process  is  evidenced  by  radial  profiles  of 
the  ambipolar  diffusion  and  net  ionization  terms  in 
the  electron  density  equation,  shown  in  Figure  9  near 
the  anode  wall  in  the  current  attachment  region. 

To  gain  an  understanding  of  the  loss  mechanisms 
involved  in  this  particular  aicjet  thruster  design,  we 
can  examine  the  partition  of  energy  m  the  flow  at 
the  exit  plane.  By  conservation  of  energy,  the  t(> 
tal  power  in  the  exiting  flow  for  the  baseline  case  is 
equal  to  the  electrical  input  power  minus  the  power 
lost  to  the  walls  plus  the  power  inherent  in  the  inlet 
flow,  for  a  total  of  14.58kW.  Table  3  catalogs  the  dis¬ 
tribution  of  this  power  in  the  various  energy  states 
of  the  exiting  flow.  As  the  table  shows,  in  this  case 


Table  3:  Distribution  of  Energy  in  the  Exit  Flow  for 
the  Baseline  Case  Arcjet  Simulation  ^  ^ 

F-nergy  State  |  Power (W7~  %  of  Total 
Axial  kinetic  5720  39.2 

Radial  kinetic _ 89 _ _ 

Azimuthal  kinetic  ~  <  1 _ 5:2 _ 

Ionization  ~  2546 _ 17^5 - 

Dissociation _ 4367 _ 30^0 _ 

Electron  thermal  ~ _ 88 _ 

Heavy  sp.  therm^  1769 _ 


only  39%  of  the  exit  plane  energy  is  in  the  form  of 
useful  thrust.  Fully  47%  of  the  energy  is  tied  up  m 
dissociation  and  ionization,  while  13%  is  classified  as 
thermal  energy.  In  particular,  the  energy  bound  in 
dissociation  and  ionization  represents  a  significant 
&ozen  loss. 


A  number  of  other  cases  have  been  run  to  eval¬ 
uate  the  model’s  effectiveness  in  predicting  arcjet 
performance  over  a  range  of  parameters.  Figure  10 
shows  simulation  predictions  of  specific  impulse  for 
the  TTl  thruster  at  three  mass  flow  rates  and  over 
a  range  of  applied  currents.  The  simulated  specific 
impulse  is  seen  to  be  approximately  5-10%  higher 
than  that  measured  by  experiment.  Table  4  com¬ 
pares  predicted  and  experimentally  measured  dis¬ 
charge  voltage  for  the  cases  presented  in  Figure  10. 
Voltage  predictions  &om  the  arcjet  simulation  fell 
from  within  1-3%  of  experiment  for  the  intermediate 
mass  flow  rate  to  within  10-15%  for  the  high  mass 
flow  rate.  The  negative  slope  of  the  V-I  relation¬ 
ship  is  also  captured.  The  differences  in  predictive 
voltage  accuracy  between  mass  flow  rates  probably 
result  from  the  boundary  condition  which  requires 
the  current  to  attach  downstream  of  the  constrictor 
exit.  In  reality  this  attachment  point  is  a  function 
of  the  operating  conditions  and  may  also  be  affected 
by  sheath  layers  or  instabilities  not  addressed  by  this 
research. 

Table  4:  Comparison  of  Discharge  Voltages  for  the 


The  effect  of  increasing  specific  power  on  arcjet 
performance  can  be  examined  in  detail  by  compar¬ 
ing  simulation  results  from  the  three  m  =  O.lg/s 
cases.  Predicted  voltages  for  these  cases  are  very 
similar,  as  are  current  patterns  in  the  nozzle.  Physi¬ 
cally,  the  distribution  of  current  is  controlled  by  the 
behavior  of  current  passage  through  the  outer  gas 
layer,  which  is  in  turn  governed  by  the  heat  diffnsiv- 
ity  in  that  layer.  Since  the  heat  diffusivity  is  essen¬ 
tially  the  same  for  these  three  cases,  the  current  pat¬ 
terns  and  voltages  are  nearly  identical  as  well.  The 
slightly  negative  V-I  characteristic  results  primarily 
from  higher  electrical  conductivity  in  the  arc  core 
due  to  increasing  electron  temperature  with  increas¬ 
ing  current  (Figure  11).  In  addition,  the  arc  width 
increases  as  the  current  is  increased  (Figure  12),  thus 
decreasing  the  width  of  the  outer  layer.  This  makes 
it  easier  for  current  to  pass  through  the  outer  layer. 


7 


$0  that  as  the  current  is  increased  the  arc  attaches 
earlier  in  the  nozzle  and  the  voltage  is  reduced. 

The  effects  of  arc  widening  and  increasing  tem¬ 
perature  with  increasing  current  result  indirectly 
from  the  greater  ohmic  dissipation  produced  by 
higher  current  levels.  This  increased  dissipation 
then  causes  additional  heating  of  the  electrons,  as 
shown  in  the  centerline  plots  of  electron  temperature 
in  Figure  11.  Higher  electron  temperatures  cause 
not  only  increased  coUisional  energy  transfer  to  the 
heavy  species,  thereby  elevating  gas  temperatures 
(Figure  13),  but  also  higher  finite-rate  dissociation 
and  ionization  rates  and  increased  ambipolar  diffu¬ 
sion,  which  cause  the  ionized  region  of  the  arc  core 
to  widen.  Interestingly,  the  electron  temperature  in 
the  flow  outside  the  main  arc  is  little  affected  by  vari¬ 
ations  in  the  applied  current,  as  shown  in  Figure  14. 
Since  the  temperature  increases  and  the  arc  core  ex¬ 
pands  radially  as  the  current  is  increased,  the  cen¬ 
tral  region  of  low  mass  flow  necessarily  deepens  and 
widens.  This  can  be  seen  in  Figure  15,  which  shows 
mass  flux  radial  profiles  at  the  constrictor  exit,  and 
in  Figure  16,  which  tracks  the  mass  flow  fractions  in¬ 
side  the  arc  as  a  function  of  axial  location.  The  inlet 
pressure  also  increases  with  the  applied  current,  ris¬ 
ing  from  l.lSatm  at  60A  to  1.37atm  at  130A.  These 
two  related  effects  produce  increased  thrust  as  the 
current  is  increased.  The  power  transferred  to  the 
anode  wall  increases  with  increasing  current,  from 
0.85kW  at  60A  to  1.98kW  at  130A,  primarily  due 
to  the  increased  flux  of  current-carrying  electrons  at 
the  anode  surface.  The  effective  anode  voltage  drop 
increases  only  slightly,  however,  from  14. IV  at  the 
lower  current  to  15.2V  at  the  higher  current. 

The  effect  of  varying  the  mass  flow  rate  on  arc- 
jet  performance  can  be  isolated  to  some  degree  by 
comparing  the  three  cases  in  Table  4  with  differ¬ 
ent  mass  flow  rates  but  similar  specific  powers.  The 
first,  fourth,  and  last  cases  listed  in  that  table  have 
specific  powers  of  120  MJ/kg  and  specific  impulses 
of  1030  seconds,  plus  or  minus  a  few  percent.  The 
main  difference  between  these  cases  lies  in  their  re¬ 
sulting  current  distributions,  as  shown  in  Figure  17, 
Since  increasing  the  mass  flow  rate  increases  the 
density,  which  directly  decreases  the  heat  diffusiv- 
ity  (a  =  — ),  the  arc  becomes  more  constricted. 
This  greater^onstriction  produces  higher  electron 
and  gas  temperatures  inside  the  arc  and  lower  tem¬ 
peratures  outside  of  the  arc  (Figure  18).  Since  the 
width  of  the  outer  layer  also  increases  with  the  mass 
flow  rate,  this  leads  to  increased  convection  in  the 
current  attachment  region  and  consequently  to  a 
broader  attachment  zone,  as  shown  in  Figure  17. 
This  effect  then  translates  into  a  higher  predicted 
discharge  voltage. 


5  Anode  Thermal  Model 

In  order  to  gain  a  more  accurate  estimate  of  the 
temperature  of  the  inlet  gas  and  of  the  anode  wall, 
an  anode  heat  balance  model  was  constructed.  Fig¬ 
ure  19  shows  a  diagram  of  the  simulated  thruster 
assembly  for  the  TTl  arejet.  The  tungsten  anode 
and  cathode,  molybdenum  casing,  boron  nitride  in¬ 
sulator,  and  propellant  gas  passages  arc  all  included 
in  this  heat  balance  model. 

The  governing  equation  in  the  thruster  assembly 
is  the  heat  flow  equation, 


where  k  is  the  local  thermal  conductivity.  For  low 
to  medium  power  arejets,  Ohmic  dissipation  is  neg¬ 
ligible  in  the  anode  but  not  insignificant  in  the  thin 
cathode  rod.  At  the  rear  of  the  thruster  assembly 
the  temperature  is  fixed  at  300®K,  while  on  the  outer 
surfaces  the  temperature  is  calculated  so  as  to  solve 
the  radiation  heat  balance 

=  eerssT^  i  (^2) 

an 

where  e  is  the  emissivity  of  the  material  and  <rsB  = 
5.67  X  10"®7^;jT  is  the  Stcfan-Boltzmann  constant. 
Along  the  inner  anode  boundary,  the  wall  heat  flux 
from  the  previously  described  flow  simulation  is  in¬ 
put,  while  at  the  cathode  tip  an  experimentally  de¬ 
termined  heat  input  of  lOOW  is  assumed. 

The  addition  of  heat  to  the  flowing  propellant  is 
essentially  a  1-D  problem,  so  the  equations 

dT  ,  8^T  '  8T  _  8^T  . 

arc  employed  in  the  gas  passages,  the  former  in  the 
axial  flow  section  and  the  latter  in  the  rewlial  flow 
section.  Parabolic  gas  temperature  and  velocity  dis¬ 
tributions  are  assumed,  and  the  equations  arc  inte¬ 
grated  by  a  4th  order  Runge-Kutta  scheme  from  the 
rear  of  the  thruster  assembly  to  the  plenum  entrance. 

For  the  case  of  (m  =  O.lg/s,  I=100A),  Figure  20 
shows  the  heat  flux  to  the  anode  as  a  function  of 
axial  distance  based  on  the  results  of  the  flow  sim¬ 
ulation  described  in  the  previous  section.  The  bulk 
of  the  heat  is  transferred  by  current  carrying  elec¬ 
trons,  which  are  deposited  in  the  first  1/3  of  the 
nozzle  in  the  current  attachment  zone.  The  resulting 
temperature  distribution  in  the  thruster  assembly  is 
pictured  in  Figure  21.  A  maximum  temperature  of 
1550®K  is  reached  in  the  current  attachment  region 
of  the  anode,  and  the  temperature  varies  from  875®K 
at  the  cathode  root  to  990®K  at  the  plenum  entrance 
to  1080®K  at  the  anode  tip.  The  propellant  in  the 
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gas  passages  is  heated  to  825'’K  by  conduction  from 
the  electrodes.  In  terms  of  power,  1.45kW  is  de- 
posited  into  the  anode  and  O.lkW  into  the  cathode. 
Of  this  input  amount,  49%  is  transferred  to  the  pro¬ 
pellant,  34%  is  conducted  to  the  backplate,  and  17% 
is  radiated  to  space.  While  beneficial  recovery  of  one 
half  of  the  anode  power  loss  is  achieved  with  this 
design,  the  results  suggest  that  performance  could 
be  increased  even  further  by  adding  thermal  insula¬ 
tion  near  the  thruster  backplate  to  maintain  a  higher 
temperature  throughout  the  thruster  assembly. 

The  resulting  inlet  gas  and  electrode  wall  temper¬ 
atures  from  the  heat  balance  model  were  then  used 
as  new  boundary  conditions  for  the  arcjet  flow  sim¬ 
ulation.  Figure  22  compares  predicted  performance 
using  both  the  old  and  new  temperature  boundary 
conditions  at  a  mass  flow  rate  of  O.lg/s.  While  the 
anode  model  shows  closer  agreement  with  experi¬ 
mental  results  at  low  currents,  it  does  worse  in  terms 
of  accuracy  at  higher  currents.  For  the  60A  case  the 
inlet  and  wall  temperatures  calculated  by  the  heat 
balance  model  are  200  -  SOO^K  less  than  those  spec¬ 
ified  in  the  original  simulation,  so  the  thrust  and  Isp 
are  lower.  For  the  lOOA  case,  the  lower  inlet  gas 
temperature  predicted  by  the  heat  balance  model  is 
offset  by  increased  heating  of  the  flow  by  the  an¬ 
ode  wall  in  the  constrictor  exit  and  current  attach¬ 
ment  regions,  resulting  in  a  slightly  higher  specific 
impulse.  For  the  130 A  case,  the  inlet  gas  tempera¬ 
ture  is  still  lower  than  originally  specified  (910®K), 
but  the  anode  wall  temperatures  are  much  higher 
{1150*K  at  the  inlet  and  1250*K  at  the  esdt,  peak¬ 
ing  at  1850®K).  This  causes  a  substantially  higher 
performance  prediction.  The  predicted  voltages  for 
these  cases  with  a  coupled  anode  thermal  model  are 
slightly  lower  than  those  calculated  with  the  previ¬ 
ously  specified  boundary  temperatures.  This  is  due 
to  preferential  concentration  of  the  anodic  arc  at¬ 
tachment  in  the  region  of  elevated  waU  temperature. 
In  summary,  use  of  the  coupled  heat  balance  model 
results  in  a  slightly  better  prediction  of  voltage  but 
no  significant  improvement  in  specific  impulse  pre¬ 
diction,  leaving  a  gap  of  5-10%  between  predicted 
and  experimental  J„.  An  anode  heat  balance  model 
such  as  the  one  described  above,  however,  does  have 
the  additional  potential  to  allow  further  optimiza¬ 
tion  of  the  thruster  design. 


6  Conclusions 

A  detailed,  multifluid,  viscous  model  has  been  d^ 
veloped  to  simulate  the  nonequilibrium  gasdynamc 
flow  in  an  electrothermal  arcjet.  The  model  is  iin- 
plemented  on  a  nonuniform  mesh  fixed  to  experi¬ 
mental  thruster  dimensions,  and  the  equations  are 


solved  using  MacCormack’s  method  and  successive 
over-relaxation.  Numerical  results  are  achieved  for 
hydrogen  propellant,  and  calculated  thrust,  spe¬ 
cific  impulse,  and  discharge  voltage  compare  well 
^th  experimental  data  for  the  German  TTl  arcjet 
thruster.  The  internal  two-dimensional  structure  of 
the  flow  is  revealed,  particularly  with  respect  to  arc 
development  and  anode  attachment,  and  the  two- 
temperature  nature  of  the  flow  is  evident.  In  par¬ 
ticular,  the  integration  of  a  separate  electron  energy 
equation  has  shown  that  stable  attachment  of  the  arc 
to  the  anode  occurs  by  increased  local  ohmic  heating 
coupled  with  nonequilibrium  dissociation  and  ioniza¬ 
tion  in  the  flow  between  the  arc  core  and  the  anode 
wall. 

An  anode  thermal  model  has  been  developed  in 
order  to  provide  more  realistic  estimates  of  the  in¬ 
let  gas  and  anode  wall  temperatures.  The  model 
predicts  a  temperature  maximum  in  the  current  at¬ 
tachment  region  due  to  local  heating  from  current- 
carrying  electrons.  Coupling  of  the  heat  balance 
model  to  the  internal  flow  simulation  results  in  closer 
agreement  with  experiment  at  low  currents  but  not 
at  higher  currents  for  a  mass  flow  rate  of  O.lg/s. 
While  the  heat  balance  model  predicts  a  lower  in¬ 
let  gas  temperature  than  originally  specified,  it  also 
shows  that  the  anode  wall  temperature  is  much 
higher  than  originally  specified  in  the  current  attach¬ 
ment  region  for  the  higher  current  cases,  resulting  in 
higher  performance  predictions.  In  summary,  use  of 
the  coupled  thermal  model  improves  voltage  predic¬ 
tions  but  overpredictions  in  voltage  and  thrust  of  up 
to  10%  still  remain.  These  overpredictions  are  most 
likely  tied  to  remaining  deficiencies  in  the  modeling 
of  the  anodic  arc  attachment. 
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Figure  1:  Electric  Potential  Axial  Profile  for  the  Baseline  Case  Aicjet  Simulation 


Figure  2:  Enclosed  Current  Contours  for  the  Baseline  Case  of  I  -  lOOA,  m  -  O.lff/s 
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Figure  3:  Ionization  Fraction  Contours  for  the  Baseline  Case  of  J  =  100>1,  m  —  O.lg/s 


Figure  4:  Electron  Temperature  Contours  for  the  Baseline  Case  of  J  =  lOOX,  m  -  O.lg/s 
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Figure  7:  Radial  Profiles  of  Some  Terms  in  the  Electron  Energy  Equation  in  the  Current  Attachment  Region 
Just  Beyond  the  Constrictor  Exit 
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Figure  8:  Radial  Profiles  of  Electron  and  Heavy  Species  Temperatures  0.25mm  Downstream  of  the  Constric¬ 
tor  Exit 
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Figure  9:  Radial  Profiles  of  Some  Terms  in  the  Electron  Density  Equation  in  the  Current  Attachment  Region 
Just  Downstream  of  the  Constrictor  Exit 


Figure  10:  Predicted  Isp  Compared  to  Experimental  Data,  German  TTl  Radiation-Cooled  Hydrogen  Arcjet 
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Figure  11:  Centerline  Ajdal  Electron  Temperature 
Profiles  for  m  =  O.lg/s  at  Three  Applied  Currents 
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Figure  12:  Ionization  Fraction  Radial  Profiles  at  the 
Constrictor  Exit  for  rh  =  O.lg/s  at  Three  Applied 
Currents 
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Figure  13:  Radial  Gas  Temperature  Profiles  at  the 
Constrictor  Exit  for  m  =  O.lg/s  at  Three  Applied 
Currents 
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Figure  14:  Anode  Wall  Axial  Electron  Temperature 
Profiles  for  m  =  O.lg/s  at  Three  Applied  Currents 
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Figure  15:  Mass  Flux  (pit,)  Radial  Profiles  at  the 
Constrictor  Exit  for  m  =  O.lg/s  at  Three  Applied 
Currents 
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Figure  16:  Mass  Flow  Fraction  Inside  the  Arc 
(a  >  0.01)  as  a  Function  of  Axial  Location  for 
m  =  O.lg/s  at  Three  Applied  Currents 
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Figure  17:  Enclosed  Current  Contours  for  £  =  120MJ/kg :  I=60A,  rh  =  O.OSg/s  (top);  I=100A,  rh  =  O.lOg/s 
(middle);  and  I=130A.  m  =  0.15g/s  (bottom) 
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Figure  18:  Radial  Profiles  of  Electron  Temperature  0.25inm  Downstream  of  the  Constrictor  Exit  for  a  Specific 
Power  of  120MJ/kg 
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Figure  19:  Arcjet  Thruster  Assembly  Schematic  for  Anode  Heat  Balance  Model 
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Figure  20:  Heat  Flux  to  Anode  from  Baseline  Case  Flow  Simulation  Results  (m  =  O.lg/s,  I=100A) 
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Figure  21:  Thruster  Assembly  Temperature  Distribution  for  Baseline  Case  Heat  Flux 


Figure  22:  Effect  of  Anode  Model  Coupling  on  Predicted  Performance  (m  =  O.lg/s) 
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A  two  dimensional  numerical  simulation  was  written  for  the  acceleration  channel  and  near  plurne 
of  a  Hall  thruster,  and  the  results  were  compared  with  experimental  measurements.  The  ^^del 
assumes  quasineutrality,  Maxwellian  electrons,  and  Bohm  diffusion  across  the  magneUc  field 
lines.  Heavy  particles  are  simulated  directly  with  a  Particle- In- Cell  (PIC)  method,  while  electrons 
are  modeled  as  a  fluid  continuum.  A  time-accurate  electron  energy  equation  is  used  to  determine 
electron  temperature,  and  a  generalized  Ohm’s  Law  is  used  to  determine  the  electric  field  strengths. 
Results  indicate  a  strong  correlation  with  experimental  performance  data.  In  particular,  the 
simulation  is  able  to  predict  thrust,  torque,  power,  and  efficiency  to  within  7%  of  the  expenmental 
values.  In  addition,  wall  erosion  rates  in  the  Russian  SPT  can  he  predicted  to  within  40%  using 
a  simple  analysis  of  the  simulation  results.  Two-dimensional  plasma  distributions  are  simdar  to 
experiment,  but  do  not  match  in  all  cases.  This  is  most  likely  due  to  the  general  Bohm  diffusion 
model  applied  globally  to  the  Hall  thruster  plasma. 


Introduction 

In  1992,  Lentz  [1]  used  a  one- dimensional  nu¬ 
merical  model  to  accurately  predict  the  oper¬ 
ating  characteristics  and  plasma  parameters  in 
the  acceleration  channel  of  a  Japanese  Type  II 
Hall  Thruster.  The  assumptions  of  the  Lentz 
model  included  quasineutrality,  Bohm  diffu¬ 
sion  across  the  magnetic  field,  constant  ratio  of 
ionization  energy  loss  to  total  electron  energy 
loss,  and  fixed  magnetic  field. 

Due  to  the  success  of  the  Lentz  one¬ 
dimensional  model,  this  research  extends  the 
model  to  two  dimensions,  using  similar  as¬ 
sumptions.  The  physical  dimensions  are  con¬ 
sidered  an  input  to  the  numerical  model,  so 
any  HaU  thruster  geometry  may  be  used,  as 
well  as  concept  designs.  A  computational  grid 
is  mapped  to  physical  space  using  nonuni¬ 
form  mapping  techniques  common  in  compu- 
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tational  fluid  mechanics.  The  magnetic  field 
is  pre-computed  from  the  iron  pole  locations 
and  solenoid  strengths.  Electrons  are  mod¬ 
eled  as  a  Maxwellian  Fluid,  while  the  heavy 
species  are  treated  with  a  modified  Particle-In- 
Cell  (modified  PIC)  methodology.  CoUisional- 
ity  is  limited  to  electron-neutral  ionization  and 
ion-neutral  momentum  exchange.  The  over¬ 
all  scheme  may  be  called  “hybrid-PIC”  since 
both  fluid  and  PIC  methods  are  used  seK- 
consistently. 

This  paper  describes  the  governing  equa¬ 
tions  and  numerical  method  used  in  modeling 
the  HaU  Thruster  in  two  dimensions.  Results 
are  then  explained,  and  comparisons  are  made 
with  experimental  data  for  Russian  SPT-class 
thrusters. 
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Governing  Equations 

0.1  Magnetic  Field 

The  magnetic  field  is  considered  equal  to 
the  vacuum  field,  and  hence  expressible  as 
B  =  Vcr.  Since  V  •  B  =  0,  it  is  also 

possible  to  define  a  magnetic  stream  function 
whose  gradient  is  everywhere  orthogonal  to  B. 
One  such  stream  function  is  A,  given  by 

=  r^  =  rBr  (1) 


dr 

dcr  p 

-7*-^  =  ^rJDz 


(2) 


(3) 


dz 

dX 

dr  ~  ""dz 
If  n  is  the  distance  normal  to  the  magnetic 
field  lines,  then, 

d  dX  d  ^d 
dn  dfi  dX  dX 

Electron  Momentum  Balance 

Along  lines  of  force,  the  electrons  are  as¬ 
sumed  to  lie  in  Boltzmann  equilibrium  with 
constant  electron  temperature 

^ -“k  ]„(„,)  =  (4) 

e 

Across  lines  of  force,  a  force  balance  on  elec¬ 
trons  may  be  written  as 

V Pe  =  -eUeE  meneUei(ui  -  Ue)  (5) 

where  Uei  is  the  electron-ion  collision  fre¬ 
quency.  Rewriting  in  terms  of  the  electron 
mobility,  fie^  and  taking  only  the  component 
normal  to  the  magnetic  field  lines. 


u 


(d4>  kT,dn,  kdTA  ,  ... 

-  =  '‘- 


(V) 


Bohm  Diffusion 

Assuming  Bohm  diffusion, 

1 

~  16B 

Applying  this  Bohm  mobility  to  Equations  (4) 
and  (6), 

d4>  _  ^  kT,dn,  rBk\Ti{n^)  dT^ 

dn  dX  eUe  dn  e  dX 

(8) 

T  d<f>‘  n  t  .  1  \ 

=  "16  IT  “  ’'7 

(9) 


Current  Conservation 

To  conserve  current,  la  =  /e  +  -ft?  indepen¬ 
dent  of  A.  Written  in  terms  of  integrals  along 
magnetic  field  lines, 

Jo  =  2xe  j  neUe^firds-2iTe  mui^ards  (10) 
Using  Equation  (8)  and  simplifying, 

dd*  _  -h  -  /o  Tie(ln(ne)  -  lyds 

dX  ^-fiSlnerHs  ^ 

(11) 


Electron  Energy 

An  electron  energy  equation  may  be  derived 
under  the  assumptions  that  electrons  have  a 
Maxwellian  velocity  distribution  and  that  the 
pressure  dyad  reduces  to  a  scalar  pressure 
term,  nekTe'. 

^  [ln,kT.'^  +  V  ■  {^n,u,kT,  +  g.)  =  5  (12) 

Here,  the  directed  kinetic  energy  of  electrons 
is  neglected,  since,  for  HaU  thrusters,  this  is 
found  to  be  smaller  than  the  random  kinetic 
energy. 

In  Equation  (12),  S  is  the  electron  en¬ 
ergy  source  due  to  ionization,  radiation,  and 
charge-field  interactions.  The  net  energy  cost 
for  producing  a  single  ion,  can  be  expressed 
as  the  sum  of  the  energy  required  for  ioniza¬ 
tion,  plus  the  energy  lost  to  excitation  of  neu¬ 
tral  atoms.  An  analytical  expression  for  (p  is 
derived  by  Dugan  et.  al.  [2].  The  result  can 
be  fitted  closely  as, 

if'  =  Ae-^  -b  C  (13) 

where  ip'  and  z  are  normalized  ion  production 
cost  and  dimensionless  electron  kinetic  tem¬ 
perature,  (p'  =  §:  and  The  constants 

A,  B,  and  C  are  given  in  Table  1.  The  vol¬ 
umetric  electron  energy  loss  rate  can  then  be 
given  as. 

Si  =  iiep'Ei  (14) 

The  electron  energy  source  due  to  the  electric 
field  is  shown  [3]  to  be, 

c  _ 

02  — 

eng 


2 


A 

B 

G 

Argon 

0.188 

0.624 

1.75 

Xenon 

0.254] 

0.677 

2.00 

Table  1:  Constants  to  fit  the  Dugan  ion  pro¬ 
duction  cost  model. 

Combining  the  sources  explicitly, 

5  =  (16) 
eUe 

Heavy  Species 

Heavy  species  are  modeled  as  discrete  par¬ 
ticles  with  negligible  temperature.  Therefore, 
only  conservation  of  mass  and  momentum  are 
applicable.  Recombination  and  charge  ex¬ 
change  are  found  to  be  small  and  both  are 
neglected.  The  ion- neutral  collision  cross  sec¬ 
tion  for  argon  and  xenon  are  assumed  to  be 
1.40  X  10“^®  and  2.15  x  10“^®,  respec¬ 
tively.  Using  these  values  with  experimental 
Hall  thruster  plasma  densities,  the  mean  free 
path  for  both  ions  and  neutrals  is  found  to  be 
large  in  most  regions.  Therefore,  ion-neutral 
momentum  exchange  is  neglected. 

The  magnetic  part  of  the  Lorentz  force  is 
ignored,  since  the  Larmor  radius  for  ions  is 
large.  Therefore,  the  force  on  an  ion  is  simply, 

Fi  =  eE  (17) 

Neutrals,  being  uncharged,  only  experience 
velocity  changes  if  they  encounter  walls.  Their 
mean  free  path  is  large  compared  to  the  scale 
of  the  device. 

Ionization  Rate 

In  this  model,  it  is  assumed  that  only 
electron-neutral  collisions  can  produce  ions. 
Also,  it  is  assumed  that  only  one  degree  of 
ionization  may  exist.  Therefore,  n*  =  TZg. 

The  nonelastic  ionization  cross  section  is 
approximated  according  to  Drawin  [4]  as  a 
function  of  electron  energy.  Assuming  a 
Maxwellian  electron  distribution,  the  resulting 
expression  may  be  written  as, 

he  =  nen„?(Te)  (18) 


where  c  is  plotted  for  argon  and  xenon  versus 
electron  temperature  in  Figure  1. 


Boundary  Conditions 

To  determine  the  electron  energy  loss  to  the 
walls,  it  is  assumed  that  the  electron  flux  is 
equal  to  the  ion  flux  normal  to  the  wall,  and 
that  all  ions  recombine  there  with  no  sec¬ 
ondary  electron  emission.  Therefore,  the  elec¬ 
tron  energy  flux  to  the  wall  is  simply, 

E^  =  2kTeTi,r,  (19) 

AH  neutrals  are  assumed  to  deflect  off  the 
walls  elastically,  at  random  angles.  Likewise, 
ions  recombine  at  the  wall  to  form  neutrals 
with  random  direction  and  velocity  magni¬ 
tudes  equal  to  the  impinging  ion’s. 

The  cathode  and  anode  are  assumed  to  be 
at  a  constant  potential  difference.  Also,  the 
electron  temperature  at  the  cathode  is  fixed, 
based  on  experimental  data,  to  2ey .  The  elec¬ 
tron  temperature  at  the  anode  is  assumed  to 
have  zero  slope.  The  inner  and  outer  walls 
of  the  thruster  are  assumed  to  be  electrically 
insulating. 

Neutrals  are  assumed  to  enter  the  acceler¬ 
ation  channel  at  the  injector.  The  injector  is 
modeled  as  an  annular  ring  in  the  center  of  the 
back  wall  (anode)  of  the  channel.  The  total 
area  of  the  injector  is  set  roughly  to  the  injec¬ 
tor  area  of  the  device  being  modeled.  How¬ 
ever,  precise  data  is  not  available  for  the  SPT, 
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so  the  injector  area  is  set  to  approximately 
30%  of  the  anode  area. 

Neutrals  are  introduced  with  an  axial  veloc¬ 
ity  equal  to  the  speed  of  sound  (y  |^)  fhe 
anode  temperature  (1000  °K  for  the  SPT). 
The  radial  velocity  component  for  each  neu- 
tral_  is  tahen  to  be  a  randoin  value  between 
-\v^\  and  juxl,  where  \v^\  =  the  mean 

directional  velocity  for  a  Maxwellian  fluid. 

The  thruster  is  assumed  to  have  an  ideal 
power  supply  which  maintains  constant  dis¬ 
charge  potential  while  allowing  the  discharge 
current  to  vary.  However,  if  the  power  level 
exceeds  a  given  point  (1  kW  for  the  SPT  op¬ 
erating  at  200  Volts),  it  goes  into  a  power  lim¬ 
iting  mode  which  drops  the  discharge  potential 
accordingly. 

Numerical  Method 

On  the  heavy  particle  time  scale,  a  time- 
accurate  solution  to  the  governing  equations 
can  be  achieved  by  iterating  successively,  as 
shown  in  Figure  2. 

A  51  X  13-node  nonuniform  spatial  grid 
is  used.  Rotational  symmetry  is  assumed. 
Therefore,  only  a  meridional  section  of  the 
Hall  thruster  is  modeled.  Grid  spacing  is  de¬ 
termined  by  the  timestep  of  the  simulation.  It 
is  set  to  the  smallest  value  which  is  much  larger 
than  the  maximum  distance  traveled  by  a  par¬ 
ticle  in  one  timestep.  Figure  3  shows  the  spa¬ 
tial  grid  used  for  the  computations  presented 
in  this  paper. 

Magnetic  Field 

The  magnetic  field  is  pre-computed  by  speci¬ 
fying  the  geometry  of  the  iron  poles  (which  are 
assumed  to  have  infinite  permeability),  and 
solving  Laplace’s  equation  on  the  regions  ex¬ 
terior  to  them.  The  method  of  red-black  or¬ 
dering  was  found  to  be  the  most  flexible.  So¬ 
lutions  take  8  hours  on  a  DEC  Alpha  worksta¬ 
tion,  but  programming  is  minimal  and  it  works 
for  a  variety  of  axisymmetric  geometries. 

Integration  of  Electron  Equations 

It  is  possible  to  combine  the  electron  energy 
equation  with  Ohm’s  Law  and  the  electron 


Figure  2:  Overview  of  the  numerical  method. 
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Figure  3:  Spatial  grid  for  the  SPT-100  geom¬ 
etry. 


current  equation.  The  result  is  of  the  form, 


dt 


+iTe  +  MTl  -VN-VOT, 


dX 
d^Te 


dX^ 


=  0 
(20) 


where  /,  J,  ...iV  represent  factors  obtained  by 
itegrating  along  magnetic  field  lines  from  the 
inner  to  the  outer  wall.  This  equation  is  a 
function  of  slow  time  parameters  and  of  the 
electron  temperature.  The  slow  time  parame¬ 
ters  are  taken  to  be  those  related  to  heavy  par¬ 
ticle  motion  {rie,  rin  and  and  those  fixed 
geometrically  (B,  r).  Holding  the  slow  param¬ 
eters  constant,  this  equation  is  solved  time- 
accurately  for  electron  temperature,  which 
evolves  faster,  using  MacCormack’s  method. 
The  space  potential  can  then  be  found  on  the 
whole  domain  by  using  Boltzmann  equihbrium 
and  Ohm’s  Law. 


Convergence 

Since  the  method  is  time-accurate,  there  is 
no  guarantee  that  the  system  will  converge  to 
a  steady  state  solution.  Two  types  of  instabil¬ 
ities  may  prevent  convergence: 

•  Plasma  instabilities  due  to  convection  and 
diffusion  of  electrons  and  ions,  combined 
with  ionization  kinetics 

•  Numerical  “noise”  due  to  the  limitations 
of  the  PIC  method  in  simulating  a  contin¬ 
uum  of  heavy  particles 

The  parameters  are  averaged  over  a  long  time 
scale,  and  iteration  is  stopped  when  the  aver¬ 
ages  reach  a  constant  value.  The  time  it  takes 
for  convergence  is  assumed  to  be,  at  the  very 
minimum,  the  time  for  convection  of  slow  neu¬ 
trals  along  the  length  of  the  grid. 

Results  and  Discussion 


PIC  Method  for  Heavy  Species 

In  this  simulation,  heavy  species  are  mod¬ 
eled  time-accurately  using  a  modified  Particle- 
In-CeU  (PIC)  method.  The  standard  PIC 
method  [5]  is  a  direct  simulation  technique 
which  models  a  gas  as  discrete  particles.  Fig¬ 
ure  4  shows  the  operation  of  one  standard  PIC 
cycle. 

The  method  used  here  has  several  unique 
features  which  enhance  performance  for  the 
specific  problem  of  simulating  a  HaU  thruster 
plasma  with  a  nonuniform  grid: 

•  Fluid  equations  and  quasineutrality  deter¬ 
mine  fields,  not  Poisson’s  equation.  This 
allows  spatial  scales  much  greater  than  the 
Debye  length. 

•  PIC  superparticles  may  have  variable 
mass. 

•  Computational  coordinates  for  each  parti¬ 
cle  are  maintained  with  a  unique  particle 
following  algorithm  based  upon  Newton’s 
method 

Complete  cases  take  20  hours  to  execute  on  a 
DEC  Alpha  workstation. 


Bishaev  and  Kim  used  diagnostic  probes  to 
measure  plasma  parameters  in  the  accelera¬ 
tion  channel  of  a  HaU  thruster.  The  results 
in  this  paper  have  been  obtained  with  a  sim¬ 
ulated  geometry  which  closely  matches  what 
they  reported  [6].  The  operating  conditions 
are  given  in  Table  2. 


Operating  Parameters  for  the  SPT 

PropeUant 

Xenon 

4>d 

200  Volts 

m 

3  mg/s 

Vacuum 

2  X  10~^  Torr 

Br^max 

.018  Tesla 

Table  2:  Operating  conditions  for  the 

Bishaev/Kim  SPT  experiments. 

The  magnetic  field  parameters  are  shown  on 
the  computational  domain  in  Figures  5  and  6. 
Although  the  radial  magnetic  field  strength  is 
close  to  the  experimental  value,  the  contours 
of  magnetic  flux  are  not  matched  exactly.  This 
is  due  to  differences  in  iron  pole  geometry,  as 
weU  as  to  the  infinite  permeabihty  and  ideal 
solenoid  assumptions  used  in  the  numerical 
model. 
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Figure  4:  PIC  computational  cycle  [5]. 
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Figure  5:  Magnetic  stream  contours,  A  [T-m^], 
computed  for  the  Bishaev/Kim  SPT  geometry. 
The  contour  numbers  correspond  to  the  values 
of  A  listed  in  the  key. 


Figure  6;  Magnitude  of  the  magnetic  field,  B 
[T],  computed  for  the  Bishaev/Kim  SPT  ge¬ 
ometry.  The  contour  numbers  correspond  to 
the  values  of  B  listed  in  the  key. 


Performance 

Although  the  Bohm  diffusion  coefficient  is 
generally  written  as, 

^  («) 

it  is  known  to  vary  by  an  order  of  mag¬ 
nitude  depending  upon  the  plasma  parame¬ 
ters.  Therefore,  the  simulation  was  run  for 
three  cases  to  isolate  the  best  fit  for  the  Hall 


thruster.  A  constant  Ksohm  was  introduced  as 
a  multiplier. 


^Bohm 


k% 

16eH 


(22) 


and  the  simulation  was  run  for  Ksohm  =  0-75, 
1.0,  and  1.25.  After  convergence,  the  results 
were  averaged  in  time.  Results  are  shown  in 
Table  3. 

It  is  interesting  to  note  that,  for  the  model 
used,  varying  KBohm  is  equivalent  to  inversely 
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varying  B.  This  is  due  to  the  form  of  the 
Bohm  conductivity.  KBohm  and  B  always  ap¬ 
pear  in  the  model  together  as  . 

The  first  case  {KBohm  —  0.75)  shows  the 
effect  of  “quenching”  the  device  with  more 
propellant  than  is  possible  to  ionize.  At  low 
KBohm,  the  power  input  to  the  electrons  is 
lower,  and  ionization  rate  is  reduced.  The 
higher  neutral  density  acts  as  an  energy  sink, 
keeping  the  electron  temperature  low  and  re¬ 
ducing  the  ionization  rate  further.  This  is  self- 
perpetuating,  since  the  ion  production  cost 
increases  at  lower  energies,  as  seen  in  Equa¬ 
tion  (13). 

In  the  second  case,  the  anode  current 
matches  the  experimental  value  almost  ex¬ 
actly.  The  efficiency  of  this  model  reproduces 
the  advertised  efficiency  of  the  SPT-100,  which 
is  50%.  The  specific  impulse  is  also  a  close 
match  with  experimental  data.  This  is  not  sur¬ 
prising,  since  the  discharge  potential  was  fixed 
to  the  same  value  measured,  and  the  simula¬ 
tion  reached  full  propellant  utilization. 

The  third  case,  KBohm  =  1-25,  represents 
the  effect  of  reduced  electron  impedance.  The 
higher  electron  mobility  increases  their  heat¬ 
ing  rate.  This  moves  the  ionization  region  to¬ 
ward  the  anode,  increasing  the  ion  wall  losses. 

The  Bohm  diffusion  mechanism  is  not  well 
understood.  It  is  an  empirical  fit  with  some 
uncertainty.  Therefore,  it  is  surprising  that 
the  standard  case  fits  the  experimental  perfor¬ 
mance  data  so  well.  Furthermore,  it  is  curi¬ 
ous  that  the  correlation  degrades  so  rapidly  as 
KBohm  deviates  from  unity,  particularly  on  the 
low  side. 

Two-Dimensional  Results 

Figures  7  through  13  show  the  time- 
averaged  two-dimensional  results  for  the  case 
KBohm  =  l-O-  Some  discrepancies  exist  be¬ 
tween  the  experimental  data  of  Bishaev  and 
Kim,  and  the  results  of  the  two-dimensional 
simulation.  The  potential  at  the  exit  of  the 
channel  is  found  experimentally  to  be  10  Volts, 
whereas  the  numerical  model  predicts  it  to  be 
118  Volts.  This  discrepancy  may  be  due  to 
the  location  of  the  cathode.  The  simulation 
cathode  is  assumed  to  be  located  about  2.3 
cm  downstream  of  the  channel.  For  the  ex¬ 


periments  of  Bishaev  and  Kim,  it  may  have 
been  located  closer.  Vacuum  chamber  back¬ 
pressure  was  thought  to  cause  this,  but  it  was 
ruled  out  by  numerically  modeling  the  neu¬ 
tral  particle  flux  from  the  chamber.  The  effect 
was  negligible,  so  chamber  back-pressure  was 
ignored.  Another  explanation  may  be  some 
type  of  enhanced  electron  conductivity  in  or 
near  the  plume  which  is  not  being  modeled. 


0.0000  0.0157  0.0314  0.0471  0.0620  0.0786  0.0043  0.1100 

2(m) 

Figure  7:  Contours  of  space  potential 

{4>  [V/m])  .  The  contour  numbers  correspond 
to  the  values  of  (f>  listed  in  the  key.  KBohm  — 
1.0;  rh  =  Z  mg/s;  <p  =  200  Volts;  K  =  3.1  A. 

The  ion  density  in  Figure  8  also  matches 
poorly  with  the  experimental  results.  The 
peak  experimental  value  is  about  7  x  10^^  m  ® 
near  the  exit  of  the  channel.  The  simulation 
predicts  it  to  be  twice  that  and  to  occur  closer 
to  the  anode. 

The  assumption  of  constant  electron  tem¬ 
perature  along  lines  of  force  seems  to  hold  only 
partially.  This  may  imply  that  the  plasma  is 
not  entirely  in  equilibrium  along  the  magnetic 
fines,  and  perhaps  the  Boltzmann  relation  is 
not  completely  applicable.  Alternatively,  per¬ 
haps  they  are  in  equilibrium,  but  an  electron 
energy  balance  holds  instead  of  Boltzmann’s 
relation.  This  is  a  matter  under  investigation. 
Nevertheless,  the  peak  electron  temperature  is 
predicted  correctly. 

Also,  the  temperature  gradients  are  strong 
in  the  experimental  measurements.  This  is 
not  the  case  with  the  simulation.  One  rea¬ 
son  for  this  discrepancy  may  be  overprediction 
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Experimental 

Numerical 

Ksohm  =  0.75  Ksohm  =  1-00 

Bohm  l*2o 

I. 

A] 

3.15 

0.39 

3.10 

2.75 

h 

A 

2.10 

0.52 

2.20 

1.90 

1st 

>  S 

1530 

200 

1475 

1400 

-F] 

'N\ 

.045 

.006 

.043 

.041 

=  ihi/m 

0.95 

0.24 

1.00 

0.86 

Va  =  hi  la 

0.67  ^ 

1.33 

0.71 

0.69 

T]e  =  {el2mi^a){F !  hy 

0.84 

0.24 

0.71  n 

0.86 

V  =  VuVaVe 

0.54 

0.08 

0.50 

0.51 

1.0 

0.03 

1.1 

2.3 

fxio  m/s] 

7.0 

1.2 

4.0 

3.1 

Torque/(F  •  r^ean) 

.02- .08 

.037 

.031 

.035 

Table  3:  Performance  results  from  the  numerical  simulation,  m  =  3  mg/s;  (f>  —  200  Volts. 
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Figure  8:  Contours  of  plasma  density 

{ui  [m”®]).  The  contour  numbers  correspond 
to  the  values  of  listed  in  the  key.  Ksohm  — 
1.0;  m  —  Z  mg/s;  <j>  =  200  Volts;  la  =  3.1  A. 


Figure  9:  Contours  of  electron  temperature 
(Te  [°-^^])-  The  contour  numbers  correspond 
to  the  values  of  Tg  hsted  in  the  key.  Ksohm  — 
1.0;  m  =  3  mg/s;  <f)  =  200  Volts;  la  =  3.1  A. 


of  the  thermal  diffusion  coefficient,  which  was 
derived  from  the  Bohm  diffusivity. 

Experimental  measurements  of  ion  flux  show 
most  of  the  flux  vectors  leaving  the  region  of 
high  plasma  density,  as  expected,  and  travel¬ 
ing  along  potential  lines  or  to  the  walls.  The 
simulation  results  (Figure  11)  show  a  similar 
pattern  of  ion  flux,  with  the  vectors  emanating 
from  a  high  density  region  closer  to  the  anode 
and  closer  to  the  inner  wall. 

The  two-dimensional  numerical  results  are 


consistent  when  taken  as  a  set.  Electron  tem¬ 
perature  increases  from  the  cathode  to  the  ac¬ 
celeration  zone,  since  there  is  little  inelastic 
collision  loss  in  that  region  due  to  the  low 
neutral  density.  Once  inside  the  channel,  the 
electrons  enter  a  region  of  higher  neutral  den¬ 
sity.  Ionization  peaks  when  the  neutral  den¬ 
sity  is  increasing.  Electron  temperature  then 
decreases  from  inelastic  losses. 
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Figure  10:  Contours  of  ionization  rate 

(nj  The  contour  numbers  corre¬ 

spond  to  the  values  of  listed  in  the  key. 
Ksohm  =  1-0;  m  =  3  mg/s;  <j)  -  200  Volts; 
la  =  3.1  A. 
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Figure  11:  Vectors  of  ion  number  flux 
(Fj  The  reference  vector  represents 

a  flux  of  14.8  X  10^“m^/s.  The  contour  num- 
bers  correspond  to  the  values  of  Vi  listed  in 
the  key.  KBohm  -  1-0;  m  =  3  mg/s;  <j)  =  200 
Volts;  la  =  3.1  A. 
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Figure  12:  Contours  of  neutral  number  density 
(n„  Ksohm  =  1-0;  m  =  3  mg/s;  (j)  = 

200  Volts;  la  =  3.1  A. 


Figure  13:  Vectors  of  ion  velocity  (ui  [m/s]). 
The  reference  vector  represents  a  velocity  mag¬ 
nitude  of  16317  m/s.  Ksohm  =  1.0;  m  =  3 
mg/s;  (j)  =  200  Volts;  la  =  3.1  A. 

Torque 

The  torque  on  the  device  may  be  written  as 
the  sum  of  the  azimuthal  moments  on  the  ions 
in  each  grid  cell  as, 

T  =  ^27rr^AceueTiiUi  x  B  (23) 

The  torque  is  generally  nonimensionalized  by 
dividing  by  the  thrust  and  the  mean  radius. 
The  torque  parameters  are  given  in  Table 
3,  and  are  consistent  with  experimental  esti¬ 
mates  [7,  8]. 
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Ion  Wall  Impingement 

From  the  numerical  simulation,  average 
wall  impingement  current  is  calculated  to  be 
1.1  Amperes.  According  to  Bishaev  and  Kim 
[6],  ion  loss  to  the  walls  is  approximately 
1.0  Ampere.  Ion  waU  currents  for  the  other 
cases  are  also  hsted  in  Table  3. 

The  SPT-100  insulator  wall  is  known  to 
sputter  most  near  the  channel  e3dt  [9].  Indeed, 
Figure  13  shows  the  mean  ion  velocities  at  the 
wall  to  be  a  maximum  near  the  exit  of  the 
channel.  At  the  inner  wall,  ion  current  density 
is  48A/m^  normal  to  the  wall.  A  comparison 
of  wall  sputtering  rate  can  be  made  by  using 
the  relation, 

6  =  kv,S,  (24) 

where  is  the  wall  reduction  rate  in  m/s, 
is  the  normal  component  of  the  ion  cur¬ 
rent  density,  and  is  the  volumetric  sputter¬ 
ing  coefficient.  Abgaryan  et.  al.  [10]  mea¬ 
sured  the  sputtering  properties  of  borosil,  a 
dielectric  wall  material  used  in  SPTs.  They 
found  that  5^  ^  .08  {mm-Y / C  for  ion  energies 
of  33  eV  at  an  80  degree  angle  of  incidence. 
Therefore, 

w  4  X  10-®m/s  (25) 

which  is  very  close  to  the  experimental  value 
of  wall  sputtering  for  SPTs.  At  the  inner 
insulator  wall,  the  initial  sputtering  rate  is 
7  X  10“®  m/s  from  experiments  by  Garner  et. 
al.  [9].  The  SPT  may  easily  have  twice  the  ion 
flux  at  that  point,  however,  in  which  case  the 
calculated  value  would  be  very  close.  Never¬ 
theless,  the  order  of  magnitude  analysis  agrees 
with  experimental  results. 

Plasma  Oscillations 

The  HaU  thruster  simulation  does  not  reach 
a  steady  state  solution.  Plasma  parameters 
fluctuate  continually,  even  after  long  conver¬ 
gence  times.  The  two-dimensional  results  pre¬ 
sented  above  have  been  averaged  over  .5  ms. 
The  anode  current  for  this  period  is  given  in 
Figure  14.  It  can  be  seen  that  fluctuates 
±15%.  The  electron  temperature  is  plotted  in 
Figure  15  for  the  same  time  period. 

Two  of  the  slow  characteristic  frequencies  in 
Table  4  can  be  seen  in  the  data.  The  low¬ 
est  frequency,  related  to  the  travel  of  neutrals 


Figure  14:  Long  time  history  of  anode  current, 
la  [Amperes]. 


Figure  15:  Long  time  history  of  electron  tem¬ 
perature,  Te  [eV],  at  z  =  1.9  cm  and  r  =  3.6 
cm. 

from  the  injector  to  the  exit  of  the  channel, 
appears  in  Figure  14  as  roughly  6  kHz.  Also, 
by  examining  Figure  15,  a  clear  28  kHz  behav¬ 
ior  is  seen,  possibly  due  to  neutrals  traversing 
the  length  of  the  ionization  region. 

To  closely  examine  the  high  frequency  os¬ 
cillations,  Ue  and  Un  have  been  plotted  on  a 
short  time  scale  (.05  ms)  in  Figures  16  and  17 
for  a  particular  node  near  the  peak  of  ioniza¬ 
tion  rate. 

The  dominant  high  frequency  information  in 
Figure  16,  may  be  seen  to  be  superimposed 
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Velocity,  m/s 

Distance,  cm 

Frequency,  kHz 

Neutrals  traverse  the  channel 

300 

4 

8 

Neutrals  traverse  the  ionization  region 

300 

1 

30 

Ions  traverse  the  channel 

3000 

4 

75 

Ions  traverse  the  ionization  region 

500  1 

1 

50 

Neutrals  traverse  two  grid  cells 

300 

.24  n 

125 

Ions  traverse  two  grid  cells 

600 

.24 

250 

Table  4:  Characteristic  frequencies  expected. 


on  a  slower  signal  which  appears  to  be  in  the 
80  kHz  range.  This  corresponds  roughly  to 
the  frequency  of  ions  traversing  the  ionization 
region. 

The  high  frequency  oscillations  visible  in 
Figure  16  occur  at  a  frequency  of  aproximately 
300  kHz.  The  neutral  density  fluctuations  in 
Figure  17  are  at  approximately  100  kHz.  Both 
of  these  high  frequency  oscillations  may  be  re¬ 
lated  to  the  time  scale  of  ion  and  neutral  pas¬ 
sage  across  grid  cells,  as  can  be  seen  in  Table  4. 

Conclusions  f 

The  Hall  thruster  simulation  accurately  pre¬ 
dicted  the  performance  parameters  for  the 
Bishaev/Kim  SPT.  Efficiency  and  specific  im¬ 
pulse  were  accurate  to  within  8%.  Anode  cur¬ 
rent  and  beam  current  were  accurate  to  within 
5%.  The  net  torque  on  the  device  was  also 
calculated,  and  was  found  to  lie  within  the 
bounds  of  experimental  measurement.  This 
implies  that  the  numerical  model  may  be  very 
useful  in  predicting  the  performance  of  alter¬ 
native  HaU  thruster  geometries. 

To  test  the  Bohm  diffusion  assumption, 
three  cases  with  different  Bohm  diffusion  con¬ 
stants  were  tried.  Surprisingly,  the  most  accu¬ 
rate  fit  was  achieved  with  a  Bohm  diffusivity 
equal  to  its  traditionally  quoted  value  of 
even  though  this  is  regarded  as  a  rough  ap¬ 
proximation  only. 

The  two-dimensional  numerical  results  fol¬ 
lowed  similar  trends  as  the  experimental  val¬ 
ues.  However,  the  simulation  predicted  higher 
ionization  rate  near  the  anode.  Also,  the  elec¬ 
tron  temperature  rise  outside  the  channel  was 
not  seen  in  experimental  results. 
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Total  wall  impingement  current  was  close  to 
the  experimental  value.  Using  the  wall  current 
density  and  an  experimental  sputtering  yield, 
the  wall  erosion  rate  was  calculated.  Compar¬ 
ing  this  to  experiments  with  the  SPT-100,  a 
close  agreement  was  found. 

Therefore,  the  numerical  model  was  useful 
for  predicting  the  performance  parameters  of 
the  Bishaev/Kim  SPT  geometry.  The  two- 
dimensional  results  lacked  close  agreement 
with  experiment,  but  were  consistent  with  the 
numerical  model.  The  usefulness  of  the  simu¬ 
lation  is  in  predicting  the  performance,  wall 
erosion  rate,  and  torque  of  prototype  HaU 
thrusters,  and  as  a  tool  for  understanding  the 
physics  of  the  plasma  acceleration  process. 

Further  Work 

Work  in  refining  the  diffusion  model  is  un¬ 
derway.  Experimental  investigations  of  the 
Hall  thruster  plasma  are  needed,  in  order  to 
generate  detailed  information  about  electron 
currents  across  and  along  magnetic  field  lines. 
This  would  allow  tabulation  of  diffusion  coeffi¬ 
cients  as  a  function  of  plasma  parameters,  and 
more  accurate  simulation  of  two-dimensional 
phenomena. 
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A  model  is  developed  to  assess  the  feasibility  of  reducing  the  frozen  losses  in  a 
hydrogen  arcjet  by  adding  very  small  amounts  of  easily  ionizable  cesium  vapor.  It  is 
found  that,  within  reasonable  constraints  on  the  constrictor  geometry,  and  without 
allowing  the  electron  temperature  to  exceed  about  7000K,  both  the  ionization  and  the 
hydrogen  dissociation  losses  can  be  essentially  eliminated,  and  a  specific  impulse  of 
about  850  seconds  can  be  obtained.  The  cesium  mass  fraction  in  the  propellant  is  about 
2%. 

Introduction. 

Frozen  losses  appear  to  be  inextricably  tied  to  the  operation  of  Electric  Propulsion 
devices:  achieving  a  sufficient  ionization  of  the  gas  already  necessitates  the  expenditure 
of  the  ionization  energy  as  a  minimum.  In  reality,  there  is  also  a  cenain  amount  of 
excitation  energy  radiated  away  (in  low  pressure  devices),  and,  of  interest  to  arcjets, 
which  tend  to  operate  with  molecular  gases,  of  dissociation  energy  as  well.  Only  a  small 
fraction  of  these  energy  investments  is  recovered  in  the  expansion  nozzle. 

An  example  from  the  modeling  work  of  Ref[l]  on  a  radiation-cooled  hydrogen 
arcjet  operating  at  .1  g/sec  flow  rate,  1.2  atm  chamber  pressure,  100  A  current  and  1 1.5 
kw  power,  shows  the  magnitude  of  the  problem:  17.5%  of  the  exit  plane  energy  flux  is  in 
the  form  of  ionization  energy,  and  30%  in  the  form  of  dissociation  energy.  Including 
thermal  energy,  this  leaves  only  39%  for  the  flow  kinetic  energy.  While  details  vary,  this 
is  typical  of  many  other  arcjets. 

A  more  subtle  form  of  inefficiency  is  associated  with  the  extreme  nonuniformity 
of  the  enthalpy  distribution  in  an  arcjet  nozzle:  the  flow  of  total  enthalpy  H  of  an  ideal 
gas  out  of  the  chamber  is  J  Hdih  while  that  of  momentum  after  full  expansion  is 
j-flHdih.  Thus,  with  no  other  losses,  the  best  possible  efficiency  is 

^  HjdmXlHdrh) 

which  is  always  less  than  unity,  the  more  so  the  less  uniformly  H  is  distributed  over  the 
flow. 
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Addition  of  a  very  small  amount  of  an  alkali  (probably  cesium)  to  the  hydrogen 
gas  in  an  arcjet  can  have  a  dramatic  influence  on  the  operating  regime  of  the  device,  and 
contribute  to  the  solution  of  both  of  the  problems  mentioned.  Since  cesium  in  a  mole 
fraction  of  a  few  times  10"^  will  be  fully  ionized  at  electron  temperatures  no  higher  than 
about  4000  K,  the  gas  can  have  an  adequate  electrical  conductivity  without  the  hydrogen 
ever  getting  ionized.  In  addition,  significant  thermal  nonequilibrium  can  develop,  so  that 
the  gas  can  be  even  cooler  than  the  electron  gas,  and  hydrogen  dissociation  can  also  be 
largely  avoided  .  This  is  limited,  of  course,  by  the  need  for  a  gas  temperature  high  enough 
for  an  attractive  specific  impulse;  however,  as  this  work  will  show,  specific  impulses  up 
to  850  sec.  can  be  easily  be  obtained  with  minimal  frozen  losses.  The  enthalpy 
distribution  is  also  much  less  nonuniform  than  in  the  classical  constricted  arc.  This  is 
because  of  the  full  ionization  of  the  seed  (cesium),  which  minimizes  the  electrothermal 
instability  responsible  for  constriction. 

The  resulting  operating  condition  bears  some  resemblance  to  that  of  a  hydrogen 
resistojet,  where  frozen  losses  are  also  minimized  by  staying  at  a  moderate  gas 
temperature.  The  difference  is  that,  while  resistojet  walls  need  to  be  hotter  than  the  gas 
itself,  the  reverse  is  true  in  the  seeded  arcjet.  Achieving  a  specific  impulse  of  850  sec  in  a 
resistojet  would  require  at  a  minimum  a  heater  wall  temperature  of  about  2300  K,  which, 
with  current  metallurgical  limits,  precludes  long  term  operation. 

We  present  in  this  paper  a  preliminary  study  of  the  performance  potential  of  the 
seeded  hydrogen  arcjet.  The  model  is  quasi  one-dimensional  and  inviscid,  and  no 
attention  has  been  paid  to  near-electrode  effects.  The  nozzle  flow  is  simply  assumed 
frozen  isentropic.  We  do  account,  however,  for  the  main  effects  in  the  constrictor,  and  are 
able  to  explore  the  practical  limits  of  the  operational  domain. 

Gas  Physics 

The  species  accounted  for  by  the  model  are  H2,  H,  Cs,  Cs^  and  e.  The  electron 
temperature  is  kept  low  so  as  to  exclude  the  presence  of  and  112"^.  Strong  coupling  is 
assumed  between  the  ions  and  the  neutrals,  designated  together  as  the  heavy  species. 
This  implies  equal  temperatures  and  velocities  (except  for  ambipolar  diffusion).  The 
electrons  are  assumed  to  be  in  Saha  equilibrium  with  cesium  atoms  and  ions 

=  (1) 

«cx  Scs  \  ) 
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where  ns  and  gs  refer  to  the  number  density  and  degeneracy  of  species  s  respectively, 
m,  is  electron  mass,  h  is  the  Planck  constant,  k  is  the  Boltzmann  constant,  and  £;  is  the 
ionization  potential  for  cesium.  The  assumption  of  Saha  equilibrium  implies  that  electron 
collisional  processes  dominate  over  radiative  processes.  This  is  true  whenever  the 
alternative  coronal  balance  model  would  predict  higher  values  of  {n^  /  which, 
using  Cs  rate  data  in  Ref[2],  would  occur  below  approxiately  /Zg  =  2x10  m  .  Our 
electron  densities  are  generally  two  orders  of  magnitude  above  this  limit.  For  the  electron 
temperatures  of  this  problem,  in  turn,  Saha's  equation  almost  invariably  predicts  full 
ionization. 

Chemical  equilibrium  is  also  assumed  to  prevail  between  hydrogen  molecules  and 
2 

atoms,  -^-K^{T)  where  K„{T)  is  the  equilibrium  constant  and  can  be  found  from 

thermodynamic  tabulations  [3].  The  electron  temperature  is  obtained  from  the  electron 
energy  equation 

cE^  =  3^6n,v,k{T,-T)  (2) 


where  cr  is  the  electrical  conductivity,  E  is  the  electic  field,  5(7 g)  is  the  inelastic 
correction  factor,  is  electron  collision  frequency  with  the  heavy  species 

The  collision  frequency  between  electrons  and  cesium  atom  is  negligible. 

The  neglect  of  radiative  losses  in  Eq  (2)  and  elsewhere  is  motivated  by  the  near 
full  ionization  of  the  seed,  since  Cs  ions  and  hydrgen  molecules  would  be  only  weakly 
excited  at  the  electron  temperatures  (under  .6  eV)  found  here. 

Equation  (2)  can  be  rearranged  as 


(3) 


where  and  depend  on  the  ratio  of  the  cross  sections  and  are  of  order  unity.  The 
above  equation  is  then  solved  iteratively  to  obtain  Te. 

Due  to  its  low  mole  fraction,  the  effect  of  the  seed  element  on  thermal 
conductivity,  enthalpy  and  density  is  negligible.  The  electrical  conductivity  is 
independent  of  pressure,  because  both,  the  electron  density  (nearly  equal  to  the  seed 
density),  and  the  collision  frequency,  are  proportional  to  pressure;  the  conductivity  is,  of 
course,  proportional  to  the  seed  fraction. 
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rioverning  Equations 

The  model  treats  a  quasi  one-dimensional  inviscid  constrictor  flow,  coupled  with 
a  radial  (Elenbaas-Heller)  diffusion  equation,  modified  to  allow  convective  energy 
transport.  The  overall  conservation  laws  are 

—  ClTtrpudr  =  0  (4) 

— (puM)  +  A^  =  0  (5) 

dx^  dx 


pu 


dh, 

dx 


r  dr\ 


'*) 


+  gE^ 


(6) 


—  llnraEdr  =  0  (7) 

dxJo 


where  A  =  nR^  is  the  constrictor  cross  section,  ht  is  the  total  enthalpy,  — 
is  the  heat  flux  potential.  Here,  Keff  is  the  effective  thermal  conductivity,  which  includes 
diffusive  transport  of  heats  of  reaction.  The  axial  electric  field  is  calculated  from  Ohm’s 
law , 


E  = 


/ 

r« 

litrcdr 

'o 


The  governing  equations  can  be  simplified  if  the  pressure,  momentum  flux,  q  =  pu^ ,  and 
electric  field  are  assumed  to  be  radially  uniform.  Expressing  all  other  functions  in  terms 
of  T  and  P,  after  algebraic  manipulation,  the  governing  equations  reduce  to 


dh,  1 
dx 


1  d  f  d(p 


rdr 


y'dr 


+  gE 


(8) 


dx 


r/?  r  fdp^ 

Jo 

dh, 

■37^'' 

ai 

^  ^  )  p  ^  ^  Jt* 

dr 

(9) 
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(10) 


dh,  dp 

yr  -T^  +  a  — 

_  =  dx  dx 

dx  p 


(11) 


where 


1 

a  = - 

2p 


^dp 


q  (dp 


Jt 


ip^ydp 


Jt 


dT 


2p^ 


ar 


These  coupled  ordinary  differential  equations  are  then  integrated  downstream  from  given 
initial  and  boundary  conditions,  until  the  flow  becomes  choked.  The  choking  point  is 
taken  to  be  the  end  of  the  constrictor,  as  well  as  the  arc  attachment  point. 

Method  of  Solution 

The  simplified  governing  equations  are  solved  for  a  fixed  constrictor 
geometry,  seed  fraction  and  current.  At  the  inlet,  the  pressure,  Mach  number  and  a 
gaussian-like  temperature  profile  with  the  peak  at  the  arc  core  are  prescribed.  An 
insulating  boundary  condition  is  specified  at  the  walls.  Other  variables  at  inlet  are  then 
computed  using  the  above  information. 


dh 

At  every  axial  station,  -z^{r)  in  equation  (8)  is  computed  as 

dx 


1 


_ 

^  ^  J.  4^ 


r.Ar^ 


(12) 


Subtituting  this  into  equation  (9)  and  integrating  along  the  radial  direction  using 

Simpson’s  method,  —  can  be  obtained.  This,  along  with  is  used  to  determine 

dx  dx 

— ,  which  then  is  used  to  determine  — .  The  step  size  is  determined  by  local  stability 

dx  dx 

and  accuracy.  Once  P,  T  and  E  have  been  determined  at  one  grid  point,  enthalpy  and 
other  variables  can  be  computed  as  well.  The  above  steps  are  then  repeated  until  sonic 
conditions  are  reached.  Arc  attachment  is  assumed  to  occur  at  the  constrictor  exit. 

Thmst  and  specific  impulse  for  a  frozen  nozzle  expansion  can  be  calculated  as 


(13) 
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where  the  thrust  coefficient,  =  1.8  (7  =  1.3  and  and  the 

stagnation  pressure  at  the  sonic  point.  The  frozen  losses  are  defined  as 


^frozen  =  A 


where  =  4AlSeV  (7. 174  x  10‘‘V)  is  the  the  dissociation  energy  of  hydrogen,  and  c  is 
the  exit  gas  velocity. 


The  model  has  been  used  for  an  exploration  of  performance  in  the  neighborhood 
of  10  kw  power.  After  some  heuristic  optimization,  a  nominal  case  was  selected,  with  the 
following  characteristics: 


Pr\  -  3.3atm; 


Molar  Seed  Fraction  =  0.0003;  1=17  A 


Constrictor  Radius=l  .4  mm;  Constrictor  inlet  Mach  No.=0.22 


The  inlet  temperature  of  the  gas  was  assumed  to  be  400K  at  the  walls,  rising  to 
3500  K  at  the  cathode  tip  (of  radius  0.35  mm)  ,  with  an  approximate  average  of  500  K. 
After  the  constrictor,  the  nozzle  was  assumed  to  expand  the  gas  to  a  Mach  number  of  5, 
with  a  ratio  of  specific  heats  of  1.3. 

The  nominal  case  was  found  to  give  these  results: 


m  =  0.369g  /  5  V=595  Volt;  Power=10.1  kw 

=  835sec;  Frozen  losses=l%;  Constrictor  L/D=8.7 

Tg  =  68(X)Ar  (Maximum) 


Computations  were  then  made  in  a  parametric  fashion,  varying  only  one  of  the 
major  parameters  at  a  time  (seed  fraction,  chamber  pressure,  current  and  inlet  Mach  no.), 
while  keeping  all  else  constant.  The  results  are  shown  in  Figures  (1)  through  (8). 

Most  of  the  trends  shown  by  these  numerical  results  can  be  understood  with 
reference  to  a  simple  model  of  heat  addition  to  a  gas  in  a  constant  area  duct.  The 
conservation  equations  for  this  situation,  written  between  the  duct  inlet  and  its  ghpked 
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exit,  are 


P*u*=PoUq^  1,  PiqMq 

]j  RgTiQ 

(16) 

P*  +  p*u*^  =  Pq  +  “  PfO 

(17) 

CpT*+^u*^  =  Cp(T,Q  +  AT,) 

(18) 

In  addition  we  must  have  u*^  =  yRgT*  and  P* /p*  -  The  approximations 
in  the  right  side  of  Eqs.  (16)  to  (18)  are  valid  for  small  Mq.  Solving  the  above  system  for 
AT; ,  we  find 

^  = 1 - ^-1  (19) 

T;0  2(r  +  l)Mo2 

which  shows  that  the  exit  total  temperature  is  proportional  to  the  inlet  temperature, 

divided  by  the  square  of  the  inlet  Mach  number  .  Since  the  thruster  specific  impulse 
scales  as  the  square  root  of  T;*,  we  find  that  I^p  should  vary  as  1/ Mq,  and  be 

independent  of  the  other  parameters.  This  is  confirmed  by  the  results  in  Figs.  ( 1 )  to  (  8  ). 
The  overall  power  balance  for  the  heated  flow  is  VI  =  mAh(  ,  which  can  be 


rewritten  as 

HA 


(20) 


With  reference  to  Eq  (19  ),  we  see  that  the  voltage  must  scale  as  the  quantity 


1 


Pro  _ 

//A'2(r  +  l)Mo 


-Mq) 


(21) 


which  is  also  confirmed  by  the  numerical  results.  In  particular,  as  Fig  (2)  illustrates,  V  is 
independent  of  seed  fraction. 

We  now  note  that  the  voltage  equals  the  mean  axial  field,  times  the  constrictor 

I  /  A 

length  L,  and,  using  Ohm's  law,  V  =  ——L.  Since  the  conductivity  is  proportional  to 


seed  fraction  and  independent  of  pressure,  we  conclude  that  the  constrictor  length  must 
scale  as  the  quantity 

Pro 


(//A)2  \y+l)MQ 


-Mq) 


(22) 


where  s  is  the  seed  fraction.  There  is  also  some  dependence  of  (7  on  the  electron 
temperature,  which  complicates  things  somewhat,  but  this  scaling  is  approximately 


validated  by  the  figures  as  well. 

Finally,  it  can  be  verified  that  the  middle  term  in  the  quadratic  Eq(  3  )  for  the 
electron  temperature  is  relatively  small,  showing  that,  to  a  first  approximation,  Tg  is 
simply  proportional  to  the  ratio  E/P  of  field  to  pressure.  Noting  that  E  is  (I  I  A)  I  g 
,  we  find  an  approximate  scaling  of  the  electron  temperature  with 
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(I  /  A)T(Q  I  {sPtoMo^).  The  dependencies  on  current  and  pressure  are  verified  by  Figs 
(4)  and  (5).  Those  on  seed  fraction  and  inlet  Mach  number  are  present  in  the  results,  but, 
due  to  the  approximations  made  in  this  discussion,  are  weaker  than  found  here. 

An  important  observation  to  be  made  is  that,  since  the  exit  gas  temperature 
depends  only  on  Mq  ,  and  the  frozen  losses  are  found  to  correspond  mostly  to 

dissociation  (which  is  governed  by  this  temperature),  these  losses  will  become  large 
below  some  critical  value  of  the  inlet  Mach  number.  Fig.  (8)  shows  that  this  happens 
below  Mq  -  0.21  for  our  conditions.  It  is  to  be  noted  also  that,  simultaneously,  the 

channel  length  is  becoming  excessive  (L/D>10,  probably  leading  to  excessive  viscous 
losses)  as  Mq  is  reduced,  and  so  is  the  electron  temperature  (  Tg  >  7000A!’  ,  probably 

initiating  ionization  of  the  molecular  hydrogen).  Increasing  pressure  can  alleviate  the 
electron  temperature  rise  (Fig  (4)),  but  it  also  aggravates  the  L/D  problem.  On  the  other 
hand,  reducing  seeed  fraction  (Figs  (1)  and  (2))  has  the  opposite  effects  of  alleviating 
the  L/D  problem  while  elevating  the  electron  temperature  further.  This  combination  of 
effects  limits  the  specific  impulse  (Fig  (8))  to  values  not  much  higher  than  that  for  the 
baseline  case. 

Some  sense  of  the  two-dimensional  degree  of  nonuniformity  of  the  gas  in  the 
constrictor  can  be  obtained  from  Figs.  (9)  and  (10),  which  show  contours  of  Tg  and 
Tg  for  the  nominal  case.  As  noted  in  the  Introduction,  there  is  much  less  arc  constriction 

than  in  a  conventional  arcjet.  Absent  from  the  calculation  are  the  wall  boundary  layers 
which  would  in  actuality  provide  a  thin  transition  to  the  wall  temperature. 

Summary  and  Conclusions 

The  results  of  the  simplified  model  presented  here  suggest  that  cesium  seeding  of 
a  hydrogen  arcjet  is  a  feasible  approach  towards  minimizing  frozen  losses  and  increasing 
the  efficiency,  while  maintaining  an  attractive  specific  impulse.  At  the  same  time,  the 
generally  reduced  temperature  levels  should  lead  to  fewer  lifetime  problems,  and  the 
presence  of  the  easily  ionizable  alkali  should  also  reduce  near-electrode  losses  (although 
these  aspects  have  not  been  analyzed). 

These  beneficial  features  are  to  be  traded  off  against  the  added  operational 
complexities  brought  about  by  the  presence  of  the  cesium.  These  include  the  cesium 
dispensing  system,  for  which  precedents  ,  including  flight  proven  designs,  already  exist, 
and  plume  contamination  concerns.  The  latter  may  be  important  in  some  applications, 
although  it  is  to  be  noted  that  a  mole  fraction  of  3  x  10“^  in  hydrogen  is  equivalent  to  a 
2%  cesium  mass  fraction  only,  which  may  begin  to  approach  the  contamination  levels 
due  to  electrode  erosion  in  other  types  of  thrusters. 
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In  view  of  these  results,  it  is  recommended  that  more  detailed  theoretical  design 
studies  be  pursued,  and  that  plans  be  made  for  experimental  work  in  this  area. 
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inlet  wall.  This  makes  the  wall  boundary  condition 
a  critical  ignition  factor.  This  problem  has  been  con¬ 
sidered  in  other  works  as  well,  many  of  which  did  not 
consider  the  details  of  the  ionising  region.  Such  mod¬ 
els  may  yield  significant  results,  but  arc  not  caF>able  of 
looking  into  failure  modes,  i.e.,  nonignition.  These  au¬ 
thors  [7],  [8]  have  previously  addressed  constant  and 
varying  temperature  models  with  constant  speed  as 
well  as  acceleration,  and  a  local  energy  balance  model 
with  accelerating  flow.  Burton  [2]  considered  a  finite 
thermal  conductivity  model,  and  included  some  radiar 
tion  effects,  with  constant  bulk  speed,  and  supersonic 
injection.  The  supersonic  inlet  assumption  essentially 
erases  any  possible  effect  that  back-diffusion  may  have 
on  inlet  ionisation,  so  that  this  should  be  considered 
as  work  of  similar  direction,  but  in  different  regimes. 

This  paper  extends  work  which  focuses  on  the  role 
that  an  upstream  flux  of  charged  particles  may  have 
on  the  gas-plasma  transition  at  the  inlet  of  a  magne- 
toplasmadynamic  thruster.  The  addition  of  an  over¬ 
all  energy  balance  consideration  to  an  accelerating 
plasma  model  makes  the  results  self-consistent. 


2  Inlet  Flow  Model 

This  section  will  present  the  basic  equations  of  mo¬ 
tion  as  well  as  details  of  the  models  used  here  for  ion 
slip  and  electrical  conductivity.  Section  3  will  consider 
ignition  in  an  accelerating  flow. 

We  will  consider  injection  of  atomic  argon  into  the 
thruster  channel  through  a  porous  backplate  to  ensure 
1-D  bulk  flow.  Figure  1  shows  the  general  configura^ 
tion  of  the  inlet  ionizing  region  (a  constant  area  chan¬ 
nel  of  width  w  out  of  the  page;  the  area  may  change 
outside  of  this  region).  Since  the  hypothesis  of  this 
work  is  that  back-diffusion  of  electron-ion  pairs  sus¬ 
tains  ignition,  the  characteristic  scale  length  of  the 
ionization  region  is  expected  to  be  the  back  diffusion 
length,  /ij,  defined  below. 
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2,1  Equations  of  Motion 

The  steady  state,  quasi-one  dimensional  equations  of 
motion  for  this  configuration  arc  overaU  and  electron 
continuity,  overall  momentum,  electron  energy,  and  a 
magnetic  field  equation.  The  equations  used  are  a 
simplified  version  of  the  1-D  equations  of  Niewood  [5]. 
The  overaU  continuity  equation  is, 


{rie  -h  na)u  = 


JlgU 


(1) 


where  Ug  is  the  total  number  density  of  of  nuclei, 
is  the  number  density  of  atoms,  and  is  the  number 
density  of  electrons.  The  mass  flow  rate  is  m,  the 
channel  area  is  ^4,  and  the  nucleus  mass  is  m,.  Since 
this  is  a  quasi-neutral  plasma,  and  only  singly  ionized 
specie  may  exist,  fie  =  n,-,  where  n,*  is  the  number 
density  of  the  ions.  The  ion  continuity  equation  is 


dneu  drieVi  Dafi^  ,  . 

H - TT-  —  *"- 

In  equation  2,  the  transverse  ambipolar  diffusion, 
is  included  (with  an  assumed  parabolic  distri¬ 
bution  of  Tie  in  the  transverse  direction  so  that  /i  is  ac¬ 
tually  of  the  actual  channel  height,  H  ^  0.02m.) 
in  order  to  allow  a  bal£uice  to  be  struck  far  down¬ 
stream  of  the  ionization  layer.  The  rate  equation  for 
production  of  ion  electron  pairs,  fie,  is 


fie  —  ^e^i'S'ea 

where  Sac  Sea  wc  the  overall  ionization  and  over¬ 
all  three  body  recombination  coefficients,  respectively. 
Values  for  these  are  given  in  references  [7]  and  [8],  and 
for  argon, 


5ea  =  8.25x  10“*^^  exp 


ToSo)  ^  3.95)^ 

0.6144 


) 
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The  integrated  total  momentum  equation,  including 
pressure  and  magnetic  (Lorenz)  forces  and  neglecting 
friction,  is 


TTiin,u^ +P+ ^  =  F  (3) 

and  the  electron  energy  equation,  including  ^aial  heat 
conduction.  Ohmic  heating,  and  energy  loss  due  to  the 
endothermic  ionization  process,  and  neglecting  both 
pressure  work  and  transverse  heat  conduction,  is 


3,  driguaTe 

r^~dr- 


(4) 


Figure  1:  1-D  Thruster  Inlet  Region 
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and  the  magnetic  field  is  governed  by  the  following 
equation 


which  is  solved  an^dytically  into  the  channel  and  to 
the  throat.  The  ionization  zone  is  found  to  be  entirely 
embedded  in  the  magnetic  diffusion  layer,  and  consists 
of  an  inner  layer,  near  the  wall,  where  ambipolar  diffu¬ 
sion  is  significant,  and  an  outer  layer,  where  diffusion 
plays  no  significant  role.  The  first  subsection  following 
describes  the  model,  and  the  second  presents  results. 
3<1  Formulation 

We  start  with  G  =  =  constant  in  the  constant 

area  ionization  zone.  The  differential  form  of  the  mo¬ 
mentum  equation  (3)  is  used: 


rriiG 


dx  dx 


P  + 


2/io , 


(12) 


where  the  pressure  for  a  two  temperature,  electrically 
neutral  plasma  is 


p  =  mingV%{Q  +  (1  -  Qr)^)  (13) 

Here,  because  of  the  acceleration,  the  alternative 
version  for  the  ion  slip  flux  is  used,  so  that,  from  cqua^ 
tion  9,  we  can  define  the  flux,  T  =  — neV^,  from 


da 

dx 


(14) 


This  equation,  along  with  the  following  three,  con¬ 
stitute  the  equations  of  motion  for  the  quasi-one- 
dimensional  isothermal  flow  problem.  Starting  with 
equation  12,  and  using  equations  13,  14,  the  momen¬ 
tum  equation  becomes  [6] 


du 

dx 


l-a(l-^)- 


Equation  2  may  be  rewritten  as: 


(15) 


and  the  final  equation  is  the  magnetic  field  equation, 
equation  5. 

Nondimensionalizing  these  equations  will  allow  us 
to  analyze  them  |>arametrically.  Defining  reference 
values  then 

the  new  variables  are  ti  =  7  =  ^  1^* 

E  =  E/{urefbo)  and  ^  The  resulting  nondi- 

mensional  equations  of  motion  are: 


-  ^i(±zlh  +  w 

_  ^  ^  IT 

■  1  -  a(l  -  0)  -  If 


5  =  ^5 


dy  da 

d?"  d? 


.  1  du 


A(l- 


(17) 

(18) 
(.19) 


db  E-xib 

d^"  l  +  (^-l)9 

where  for  convenience  we  define 

,  i  +  (i-i)9_  (dby^ 

E  —  ub 


(20) 


In  addition  to  the  ]>arameters  defined  previously, 


{Ca)f^o^ref^ref  ^  (^<i)|| 

G  ^  Am 


A  = 


SqjG 


_  CqAm 

The  first  far  downstream  (^  oo)  boundary  con¬ 
dition  will  be  that  E  =  u6,  which  implies  an  infinite 
length  and  hence  an  infinite  magnetic  Reynolds  num¬ 
ber.  This  motivates  us  to  choose  the  magnetic  field 
itself  as  the  independent  variable.  Dividing  the  three 
other  equations  of  motion  (17,  18,  19)  by  equation  20, 
we  arrive  at  the  following  set  of  equations: 

du  mi-e):^-2eb 
db  -  i_a(l-d)-i4 


The  second  downstream  (oc)  boundary  condition  is 
the  same  balance  between  transverse  ambipolar  diffu¬ 
sion  and  ionisation  used  in  the  constant  speed  case, 
which  is  now: 

^  SionG^h^  ^  A 
Additional  boundary  conditions  arc  the  conditions 
at  the  injector  wall,  and  the  internal  boundary  con¬ 
ditions  at  the  sonic  passage  point.  The  injector  wall 
condition  is  critical:  the  ions  (moving  at  a  velocity 
equal  to  the  average  velocity  plxis  the  ion  slip  veloc¬ 
ity)  must  reach  the  sheath  at  the  inlet  wall  at  the 
Bohm  velocity  moving  upstream  against  the  bulk  flow: 


nc(ti  -h  V;)  =  -ri^VB 
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approaches  that  of  the  magnetic  scale  «  1),  At 
higher  A,  the  ionising  scale  is  decidedly  smaller  than 
the  magnetic  scale,  approaching  ^  «  fx?.  As  found 
with  the  constant  speed  case,  again  higher  A  yields  a 
higher  a^. 


therefore  makes  it  self-consistent  and  will  indicate  a 
more  definite  ignition  criterion. 

Assuming  infinite  thermal  conductivity,  which  im¬ 
plies  a  constant  electron  temperature,  and  is  consis¬ 
tent  with  the  constant  temperature  assumption  made 
previously,  the  integrated  energy  equation,  as  reported 
by  Lawless  and  Subramanian  [3]  is 
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Figure  3:  Ionisation  fraction  a  vs  ^  for  the  “standard” 
parameters,  except  for  varying  A. 

Using  equation  13  in  the  integrated  momentum 
equation  (equation  3),  and  with  constant  G  = 
in  the  ionisation  sone,  and  using  the  nondimensional 
variables  and  parameters  defined  above, 


which  is  valid  for  the  constant  area  region  assumed 
for  the  ionisation-acceleration  sone.  C  is  a  constant 
of  the  flow,  and  the  enthalpy,  h,  is  defined  as 

\2  mi  mi/ 

for  fcoien  heavy  species  temperature,  T,.  The  en¬ 
ergy  balance,  for  constant  temperature,  and  using  the 
nondimensional  variables  and  parameters,  is  now: 

2£(l-6oo)  =  (ttoo-Of.)  + 

where 


.  /?(tt  -Kl  -  a)9)  ^  ^  ^  (29) 

ti  miCUrtf 

For  small  (“subsonic”)  speeds,  <  p{oc  +  (1  -  a)^)» 
the  acceleration  comes  mostly  from  the  pressure  gra¬ 
dient  (which  is  strongly  dependent  on  the  gradient  of 
a  in  this  model),  so  that 

g(»  +  (l-»w..  f  (31,, 

ti  TTliGUreJ 

which  depends  strongly  on  the  value  of  A.  For  large 
ti  (“supersonic”),  the  acceleration  comes  chiefly  from 
the  changing  magnetic  field: 

a  +  (31) 

m»GlXrc/ 

which  does  not  depend  very  strongly  on  A  so  long  as 
the  entire  ionisation  region  is  Coulomb  dominated. 

3.3  Overall  Energy  Balance 


In  reality,  the  electron  temperature,  and  therefore  the 
parameters  of  the  problem  (namely  A,  will  be  func¬ 
tions  of  the  expansion  ratio.  This  means  that,  for 
fixed  channel  throat  ^mensions  and  flow  rate,  there 
may  be  only  one  possible  A  for  each  expansion  ratio. 
This  is  required  by  the  energy  balance  in  the  ioniza^ 
tion  sone.  Adding  an  energy  balance  to  this  model 


and  5  =  0.59  for  the  argon  atom,  and  the  Ure/  defined 
here.  This  depends  on  conditions,  of  course.  In  fact, 
{  1  is  a  rough  onset  criteria;  i.e.,  the  flow  kinetic 

energy  approaches  the  ionisation  potential  energy  of 
the  atom. 

If  aoo>ao,  and  Uoo>tio»  which  has  been  generally 
found  in  the  results  of  the  accelerating  ignition  model, 
then 

2£(1  -  6«,)  =  aoo  (1^(1 -<>)  +  «)  + ^ 

Since,  from  the  far-downstream  condition, 
a«  =  l-  — 

and  tioo  =  E/boo 


or,  rearranging. 


Which  i*  now  a  function  of  the  temperature- 
dependent  j>arametert,  and  the  (£,  600)  eolution  pair. 
Thi*  condition  sets  the  electron  temperature,  and  an 
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Abstract 

An  overall  ionization  rate  model  has  been  de¬ 
veloped  which  takes  into  account  the  detailed 
effects  of  processes  involving  excited  states, 
and  which  may  be  used  in  a  variety  of  flow  mod¬ 
els.  The  resulting  overall  collisional  rate  coef¬ 
ficients  for  3-body  recombination  are  reported 
here  for  the  argon  atom  and  first  ion,  and  for 
the  hydrogen  atom.  The  phenomenon  of  inlet 
ionization^  ignition  (steady  state  initiation  of 
ionization)  in  self-field  MPD  thrusters  has  also 
been  further  analyzed  and  explained.  The  hy¬ 
pothesis  that  back  diffusion  could  explain  the 
mm-scale  ionizing  regions  found  in  experiments 
has  been  tested  in  constant  speed  models  of 
atomic  injection  (first  without  and  then  with 
variable  temperature)  along  with  an  accelerat¬ 
ing,  subsonic  atomic  injection,  constant  tem¬ 
perature  model,  and  can  explain  ignition  over 
a  range  of  conditions  identified  here. 

Nomenclature 

A  channel  area  (A  =  Fn;). 

B  magnetic  field, 

c  a  thermal  speed. 

Ca  modified  ambipolar  coefficient  =  D^ng 
Da  ambipolar  diffusion  coefficient, 

e  electron  charge. 

E  electric  field. 

Eat  ionization  potential  energy. 

G  flow  rate  per  unit  area  =  n^t*. 

H  channel  height. 

h  channel  height/ >/l2. 

ks  Boltzmann  constant. 

I  a  scale  length. 

Kt  electron  thermal  conductivity. 
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L  thruster  characteristic  length. 

Le  Lewis  number. 

m  mass  flow  rate  =  miTiguA. 

m  particle  mass, 

n  particle  number  density. 

p  pressure. 

q  critical  value  of  a  in  conductivity  equation. 

Rfn  magnetic  Reynolds  number,  Rm  =  p><tuL, 
Sjk  collisional  rate  coefficient,  j  ^  k  transition. 
T  particle  temperature  (K). 

It  plasma  bulk  velocity. 

V  slip  speed. 

z  streamwise  coordinate. 

a  ionization  fraction,  a  = 

0  =  vl/Un. 

r  ion  slip  flux:  T  =  —Tie  Vi 

6  the  temperature  ratio: 

Am  magnetic  interaction  length. 

<r  scalar  conductivity, 

r  nondimensional  electron  temperature. 

<f>  potential  (voltage). 

SubscripU/superscripU: 

D  diffusion, 

c  electrons. 

g  heavy  gas  species. 

i  ions. 

o  inlet. 

$  sonic  passage. 

i  throat. 

.00  far  downstream  asymptotic  value 
♦  equilibrium  value. 


1  Introduction 

The  fust  section  of  this  papei  presents  the  results 
from  a  method  of  calculating  overall  recombination 
and  ionization  coefficients  for  both  atomic  and  ionic 
species.  Results  are  given  for  argon  atoms  and  ions 
and  for  the  hydrogen  atom. 
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The  rest  of  this  paper  focuses  on  the  role  that  an 
upstream  flux  of  charged  particles  may  have  on  the 
gas-plasma  transition  at  the  inlet  of  a  magnetoplasma- 
dynamic  thruster.  Steady-state  “ignition”  is  defined 
here  as  the  condition  when  this  ionizing  zone  is  fixed 
at  the  inlet  waU,  making  the  wall  boundairy  condition 
a  critical  ignition  factor. 

This  problem  has  been  considered  in  other  works 
as  well,  although  the  details  of  the  ionizing  zone  are 
generally  ignored.  Such  models  do  yield  significant  re¬ 
sults,  but  are  not  capable  of  looking  into  failure  modes, 
i.e.,  nonignition.  These  authors  [7],  [9]  have  previously 
addressed  constant  temperature  models  with  constant 
speed  as  well  as  acceleration,  and  a  local  energy  bal¬ 
ance  model  with  accelerating  flow.  Burton  [2]  consid¬ 
ered  a  finite  thermal  conductivity  model,  and  included 
some  radiation  effects,  with  constant  bulk  speed,  and 
supersonic  injection.  The  supersonic  inlet  assump¬ 
tion  essentially  erases  any  possible  effect  that  back- 
diffusion  may  have  on  inlet  ionization,  so  that  this 
should  be  considered  as  work  of  similar  direction,  but 
in  different  regimes. 

Three  cases  will  be  addressed  here:  atomic  injec¬ 
tion  at  constant  speed  with  constant  temperature 
and  then  varying  temperature,  and  atomic  injection 
with  constant  temperature  and  accelerating  flow.  All 
units  herein  are  mks,  and  temperatures  are  in  degrees 
Kelvin. 


at  an  overall  recombination  coefficient  (useful  up  to 

T,  =  500Q0K): 


8,25x10“^^  exp 


(' 


0.6144 


(1) 


This  compares  favorably  with  the  results  reported 
by  Owano,  et  al  [6],  which  adjusted  calculated  over¬ 
all  rate  coefficients  (based  on  experimentally  derived 
cross-sections)  with  spectroscopically  measured  rates 
in  argon,  and  arrived  at  the  following  overall  three- 
body  electron-ion  recombination  coefficient: 


sfr"’  =  3.ax  i0-«  +  2)  .IP  ^ 

(2) 
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Owano 

6.10  X  10-« 

3.15  X  10-« 

5.87  X  lO-*^ 

Table  1:  Collisional  rate  coefficients  for  three-body 
recombination  in  the  argon  atom  (AI):  comparison  of 
the  results  of  this  work  and  those  of  Owano.  Units  of 
Sea  m^/s. 


2  The  Ionization  Rate  Model 

This  section  extends  the  work  reported  by  these  au¬ 
thors  previously  on  the  subject  of  calculating  overall 
ionization  eind  recombination  coefficients  for  atomic 
and  ionic  species.  The  formulation  has  been  presented 
[8]  and  here  results  are  given  for  the  argon  atom  and 
ion  and  the  hydrogen  atom. 

In  most  of  the  channel,  the  excited  states,  which 
are  relatively  close  energeticafly,  may  be  expected  to 
relax  fast  enough  for  reactive  balance  to  be  assumed. 
This  dynamic  equilibrium  of  the  excited  states  holds, 
regardless  of  what  the  plasma  is  doing  if  the  excited 
state  time  scales  are  short  enough.  Given  the  excited 
state  population  distribution,  and  assuming  that  the 
ground  state  population  is  essentially  the  neutral  den¬ 
sity,  then  overall  rate  coeflficients  may  be  calculated. 
Presented  here  are  overall  recombination  coefficients. 
The  overall  collisional  ionization  coefficient  may  be 
calculated  by  imposing  overall  microreversibility  [8]. 

2.1  Argon  Atom  (AI) 

The  argon  atom  model  consists  of  21  levels,  including 
the  ground  state,  19  lumped  excited  states,  and  the 
continuum  (the  ion).  The  results  using  the  Drawin 
cross-section  model  [4]  have  been  curve-fitted  to  arrive 


2.2  Argon  Ion  (AH) 

A  33-level  second  ionization  model  for  argon  was  also 
developed.  Note  that  the  gap  between  the  highest 
excited  state  used  and  the  ion  is  2.77  eV.  This  is  a 
weeikness  at  low  Te  (as  the  low  temperature  behavior 
depends  on  the  missing  upper  levels)  and  will  have  to 
be  accounted  for. 

3485)^  j 

(3) 

2.3  Hydrogen  Atom  (H) 

The  hydrogen  atom  has  been  modeled  with  19  levels. 
The  resulting  collisional  recombination  rate  coefficient 
is  (good  up  to  Te  =  60000K): 

S„  =  6.985  X  10-«exp  ( MMzilHHl!)  !!!! 

\  0.8179  )  $ 

(4) 


Sec  =  7.17  X  10 
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exp  — 
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0.9293 
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3  Inlet  Flow  Model 


The  rest  of  this  paper  will  investigate  steady-state 
ionization  at  the  inlet,  or  ignition,  of  a  self-field 
magnet oplasmadynamic  thruster.  This  section  will 
present  the  basic  equations  of  motion  as  well  as  de¬ 
tails  of  the  models  used  here  for  ion  slip  and  electrical 
conductivity.  Section  4  will  consider  three  cases  of 
constant  speed  ignition,  and  section  5  will  consider 
ignition  in  an  accelerating  flow. 

We  will  consider  injection  of  atomic  argon  into  the 
thruster  channel  through  a  porous  backplate  to  ensure 
1-D  bulk  flow.  Figure  1  shows  the  general  configura¬ 
tion  of  the  inlet  region.  Since  the  hypothesis  of  this 
work  is  that  back-diffusion  of  electron-ion  pairs  sus¬ 
tains  ignition,  the  characteristic  scale  length  of  the 
ionization  region  is  expected  to  be  the  back  diffusion 
length,  /u,  defined  below. 
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as  ^  ♦  oo 
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Figure  1:  1-D  Thruster  Inlet  Region 

3.1  Equations  of  Motion 

The  steady  state,  quasi-one  dimensional  equations  of 
motion  for  this  configuration  are  overall  and  electron 
continuity,  overall  momentum,  electron  energy,  and  a 
magnetic  field  equation.  The  equations  used  are  a 
simplified  version  of  the  1-D  equations  of  Niewood  [5]. 
The  overall  continuity  equation  is, 

{tie  +  na)u  =  UgU  =  G  =  (5) 

Anti 

where  Ug  is  the  total  number  density  of  of  nuclei, 
is  the  number  density  of  atoms,  and  is  the  number 
density  of  electrons.  The  mass  flow  rate  is  m,  the 
channel  area  is  A,  and  the  nucleus  mass  is  n^.  Since 
this  is  a  quasi-neutral  plasma,  and  only  singly  ionized 
specie  may  exist,  n*  =  n^,  where  rii  is  the  number 
density  of  the  ions.  The  ion  continuity  equation  is 

drteU  dn^Vi  . 

-sr=—z, — S5-+".  («) 

In  equation  6,  the  transverse  ambipolar  diffusion, 
is  included  (with  an  assumed  parabolic  distribu¬ 


tion  of  Tie  ^  transverse  direction  so  that  h  is  actu¬ 
ally  of  the  actual  channel  height,  H  «  0.02m.)  in 

order  to  allow  a  balance  to  be  struck  far  downstream 
of  the  ionization  layer. 

The  integrated  total  momentum  equation,  including 
pressure  and  magnetic  (Lorenz)  forces  and  neglecting 
friction,  is 


TTljTljU* +P+  ^  f  (7) 

and  the  electron  energy  equation,  including  axial  heat 
conduction,  Ohmic  heating,  and  energy  loss  due  to  the 
endothermic  ionization  process,  and  neglecting  both 
pressure  work  and  transverse  heat  conduction,  is 


3,  dugUaTt  _  d  f  ^  dTe\ 
2*^  dx  “  dx  dx  ) 


(8) 


and  the  magnetic  field  is  governed  by  the  following 
equation 

jp 

—  =  -fXoj  =  -fio<r[E  -  uB]  (9) 

where  is  the  electron  (and  ion)  number  density  and 
Tin  is  the  neutral  particle  number  density,  u  is  the  bulk 
speed,  p  is  the  pressure,  B  is  the  magnetic  field,  T*  is 
the  electron  temperature,  E  and  j  are  the  electric  field 
and  current  density,  and  G  and  F  are  constants  of  the 
flow.  Vj  =  Uj  —  ti  is  the  species  j  slip  speed,  and  Uj 
is  the  speed  of  species  j.  The  ionization  fraction  is 
defined  as  a  =  n^jug. 


3.2  Ion-Neutral  Slip 

Since  ions  are  being  directly  accelerated  while  neu¬ 
trals  are  accelerated  by  collisions  with  ions,  there  are 
actually  two  effects  (acceleration  and  density  gradi¬ 
ents)  leading  to  ion-neutral  slip.  In  [9]  we  derived  an 
alternative  to  Fick’s  Law  from  the  neutral  momentum 
equation.  Using  the  definition  of  the  ambipolar  diffu¬ 
sion  coefficient: 


D,  =  !isi2i+l!lA 


as  well  as  the  definition  of  the  Bohm  velocity: 


^  W  +  T,) 

V  ^ 


and  the  temperature  ratio 


(10) 


(11) 


(12) 


T. +  T, 

[6  should  be  small  near  the  inlet,  and  we  use  ^  =  0.1 
here)  we  arrive  at  the  following  expression  (from  the 
neutral  momentum  equation)  for  the  ion  slip  flux: 
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n.v;-  =  D,n,e  (il  -  a) 


1  du 
u  dz 


dz) 

(13) 


If  there  is  no  acceleration  and  no  total  pressure  gra- 
dient  {0^^  =  “’^)j  Fick’s  law  results: 


dn»  dot 

{n,V,)jr  =  =  -Dan,—  (14) 

The  ambipolar  diffusion  coefficient  can  be  replaced 
by  Ca  =  DaUg^  a  function  of  temperature  only.  Using 
the  form  for  Da  iu  equation  10,  Ca  is,  for  argon 


(7.  =  1.274*10*® 


3.3  Electrical  Conductivity 

The  electrical  conductivity  which  appears  in  equations 
8  and  9  is  significant  for  both  the  Ohmic  heating  and 
the  rate  of  change  of  the  magnetic  field,  and  is,  in 
general,  a  function  of  both  electron  temperature  and 
ionization  fraction.  This  behavior  is  seen  at  low  a, 
and  may  be  expressed  as  [1]: 


-f  I^en) 

where  the  electron-ion  collisional  rate  goes  as  the  elec¬ 
tron  temperature  to  the  —3/2  power: 

aTr®/*ln  ^1.24  X  10*' 

using  mks  units,  and  temperatures  in  K,  and  the 
electron-neutral  collisional  rate  goes  as  the  square  root 
of  the  electron  temperature  [1]: 


q  is  roughly  the  minimum  ionization  fraction  for 
Coulomb  collision  dominance  at  the  reference  temper¬ 
ature.  q  0.001  results  from  typical  channel  con¬ 
ditions.  The  reference  conductivity  is  approximately 
1000  to  4000  Si. 


4  Constant  Speed  Ignition 

In  this  case,  we  extend  the  work  in  [9],  adding  the 
effects  of  temperature  variation  to  a  constant  speed 
ionizational  ignition  analysis  in  a  constant  area  chan¬ 
nel.  The  total  number  of  nuclei  (which  is  a  constant 
for  constant  speed)  in  the  plasma  is  governed  by  equa¬ 
tion  5,  with  G  constant. 

Using  Tick’s  law  in  this  case,  equation  14,  then  the 
definition  of  the  ion  slip,  T  is: 

da  ^  UeVi  _  r 
dx  Ca  Ca 

Taking  advantage  of  the  constant  G,  the  convection 
term  in  equation  6  becomes  =  GT/Ca^  by  equa¬ 
tion  16,  so  that  the  ion  continuity  equation  is  now 

dr  Gr  , .  Dane  ^ 

*  =  c:  -  <”> 

Equation  8,  the  energy  equation,  is  modified  by  the 
assumption  that  j  «  and  defining  the  heat  flux 
as 


q^Ke 


dx 


(18) 


so  that  the  lesnlting  energy  equation  is 


3  doT. 


dx 


-j-  —  Eaine  (^^) 


The  nondimensional  forms  of  equations  16,  17,  18, 
and  19  are,  respectively: 


^en 


oc  r*/® 


Assuming  that  In 


10,  and  using 


numerical  values  for  argon  from  Bittencourt  [1],  then 
the  electrical  conductivity,  nondimensionalized  by  a 
reference  value,  is 


_ ,15, 

a  +  (1  -  a)qT^/^  1  +  (^  -  1)9T*/* 

where  t  =  Tt/Trtf  and  T,,/  is  the  electron  tempera¬ 
ture  as  z/lo  —*  CO,  and 

g  =  7.2  X 

=  7.76  X  Si 


da 

di 


=  7 


(20) 


=  7  -  aAi  [5(1  -  a)  -  a* A,  -  A3]  (21) 

%  =  Q  (22) 

^  =  -^[aQ  +  r7  -  o-n  -f-  ©«jaAi[5(l  -  a)  -  a^Ajj] 
of  Le 

(23) 

where  the  length  is  nondimensionalized  as  ^  =  x/lj)^ 
where  Ij)  =  Ca/G  is  the  back-diffusion  length,  the 
nondimensional  slip  is  7  =  F/G,  and  the  electrical 
conductivity  is  ^  =  aftTref^  The  nondimensional  pa¬ 
rameters  are: 
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extinction 


Ca<rr,fE^ 

^kBT„fG^ 


Ai  = 


GaS^i 

~iir~ 


As  = 
©«  = 


Sai{hngy 

E^i 

^ksTref 


The  Damkohler  coefficient  for  ionization  (Ai  = 
hlhi)^  The  nondimensional  ionization  coefficient  is 
S(t)  =  and  is  evaluated  at  Tref,  the 

electron  temperature  at  ^  =  oo.  The  Lewis  num¬ 
ber,  which  is  the  ratio  of  the  diffusion  length  scale  to 
the  thermal  conductivity  length  scale  and  is  typically 
O(10~^),  is 

3  k^rigDa  _  Ijy 
"^2  Ke 

4.1  Atomic  Injection,  Constant 

The  case  of  constant  speed  and  temperature,  and 
atomic  injection  was  treated  previously  by  these  au¬ 
thors  [9],  and  is  repeated  briefly  here  since  its  resTilts 
introduce  some  points  relevant  to  the  other  cases.  Set¬ 
ting  (3  =  0,  leaves  only  two  equations,  for  a  and  y 
(equations  20  and  21).  Using  the  ionization  fraction 
as  the  independent  variable  (dividing  equation  21  by 
equation  20,  the  result  is 

il  =  1  _  Aia[(l  -  g)  -  a^Aj  -  A3] 
da  y  ^ 

The  first  boundary  condition  chosen  is  critical;  that 
the  ions  (moving  at  a  velocity  equal  to  the  average 
velocity  plus  the  ion  slip  velocity)  reach  the  sheath  at 
the  inlet  wall  at  the  Bohm  velocity  moving  upstream 
against  the  bulk  flow: 

n.(tt  +  Vi)  = -n.t»s 

Using  the  nondimensional  variables,  this  becomes 

To  =  ao  (l  +  (25) 

The  second  boundary  condition  is  that  the  deriva¬ 
tive  be  smooth  as  7  0,  so  that  the  numerator  of  the 

fraction  in  equation  24  is  zero  at  ^  =  00: 

1  —  ttoc  —  —  As  =  0  (26) 

so  that  the  ionization  faction  at  00  is 

a<x>  =  [v/l  +  Aj(l  -  As)  -  ij  (27) 

The  ignition  criterion  is  found  by  a  local  analysis 
of  equation  24  near  a  =  7  =  0.  Successful  ignition  is 
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Figure  2:  Blowoff  speed  for  argon,  in  y,  constant 
speed  and  temperature,  atomic  injection. 


defined  as  as  a  solution  where  the  wall  condition  may 
be  satisfied  at  7©  >  0  for  a  finite  ionization  zone  [7], 
requiring 

Ai(l-A3)>i  (28) 

which  sets  a  lower  limit  on  the  Damkohler  coefficient, 
or  an  upper  limit  on  the  injection  speed.  This  upper 
limit  speed  can  be  considered  a  ‘‘blowofP  speed; 


.  <  2,  -^^  =  ...  (29) 

This  blowoff  speed  is  plotted  vs  Te  in  figure  2  for 
various  G  and  JST  =  0.02m.  In  this  case,  we  use  the 
argon  ionization  coefficient  described  earlier  in  this  pa¬ 
per  which  differs  from  figure  1  of  reference  [9],  which 
has  a  blowoff  speed  using  the  Hinnov-Hirschberg  ion¬ 
ization  coefficient. 

Therefore,  successful  ignition  in  this  simple  model 
sets  an  upper  bound  on  the  injection  speed,  as  is  found 
in  simple  1-D  diffusion  flames.  Note  also  that  G  ap¬ 
pears  in  the  product  A1A3  in  equation  28,  so  that, 
under  certain  conditions,  where  the  00  boundary  con¬ 
dition  is  accepted,  there  is  also  a  constraint  of  mass 
flow  rate,  m  =  m<GA.  It  should  be  noted  that  in  the 
cases  displayed  above,  the  blowoff  speed  is  less  than 
the  acoustic  speed  of  the  heavy  particles, 

4.2  Atomic  Injection^  Varying 

In  this  section  we  now  allow  T*  to  vary  (finite  thermal 
conductivity)  to  investigate  the  effects  on  the  ignition 
criterion.  Using  a  as  the  independent  variable  again, 
the  equations  for  constant  speed,  finite  thermal  con¬ 
ductivity,  and  atomic  injection  are  equations  21,  22, 
and  23,  respectively,  all  divided  by  equation  20: 
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(31) 


da 


1^ 

Le 


T  + 


aQ  “  ^-n  +  ©aiaAi[5(l  -  a)  ~  a^A2] 


(32) 

The  boundary  conditions  are  again  the  wall  condi¬ 
tion,  equation  25,  plus  that  aU  derivatives  are  smooth 
at  the  asymptotic  limit:  7oo  0.  One  such  condition 
will  be  the  same  as  equation  26.  In  addition,  Qoo  =  Oj 
and 


Figure  3:  Constant  speed  ignition  with  atomic  injec¬ 
tion  and  varying  temperature:  a  and  T®. 


cE  —  ©ataAi[l  —  aoo  —  =  0  (33) 

Finally,  it,  is  assumed  that  the  electrons  do  not  heat 
the  inlet  wall,  so  that  the  initial  electron  heat  flux  is 
zero:  Qo  =  0. 

Note  that  5  =  5'=lat^  =  oo,  since  that  is  where 
the  reference  values  are  taken.  The  combination  of 
equations  26  and  33  is  usually  used  to  solve  for  aoo 
and  =  T^tf  as  a  function  of  (n^,  h,  Tj,  ^?).  In 
this  case,  the  specified  values  are  (G,  Trc/jUi  A)  {u  is 
arbitrary  since  no  momentum  balance  is  used),  so  that 
(aooi  E)  =  f{ng{u,  G),  h, Zef),  and 


n  =  ©ctAiAattoo 

Local  analysis  near  a  =  7  =  0  results  in  the  a  smulaj 
ignition  criterion  as  in  the  constant  Te  case,  except  it 
is  evaluated  at  the  inlet  temperature: 


E  = 


EaiCaao 


^ref 


u<2 


CgSaijO)  _ 


=  ttfco 


(34) 


Figure  3  shows  the  temperature  and  ionization  pro¬ 
files  for  the  trajectory  corresponding  to  Le  =  0.037, 
Ai  =  0.285,  n  =  0.56.  The  corresponding  dimensional 
quantities  are  denoted  on  the  graph. 

It  appears  that  the  temperature  variation  does  not 
have  a  dramatic  qualitative  effect  on  ignition.  How¬ 
ever,  it  will  have  a  quantitative  effect,  particularly 
at  temperatures  below  lOOOOif,  where  the  ioniza¬ 
tion  rate  coefficients  are  very  sensitive  to  tempera¬ 
ture  changes.  Another  way  that  the  temperature  effect 
may  be  greater  would  be  when  the  Hall  effect  plays  a 
role  in  the  thermal  conductivity,  raising  the  inlet  tem¬ 
perature  even  higher,  while  shrinking  the  zone  where 
temperature  variations  are  significant.  Burton  [2]  con¬ 
sidered  this  effect.  Note  also  that  Te  at  the  wall  is  al¬ 
ways  greater  than  T«  at  00,  so  that  the  blowoff  speed 


in  this  case  wiU  always  be  higher  than  that  for  the 
constant  temperature  case. 


5  Ignition  in  Accelerating 
Flows 


We  now  drop  the  assumption  of  constant  speed, 
with  Q  =  0.  This  allows  us  to  analyze  the  interac¬ 
tion  between  the  behavior  of  the  magnetic  field  and 
the  ionization  process.  The  problem  is  treated  as  two 
inner-outer  problems.  First,  there  is  a  magnetic  inner 
(diffusive)  region  &om  the  wall  out  to  a  zero  current 
boundary  condition,  and  then  an  ideal  outer  region 
which  is  solved  analytically  into  the  channel  and  to 
the  throat.  The  ionization  zone  is  found  to  be  entirely 
embedded  in  the  magnetic  diffusion  layer,  and  consists 
of  an  inner  layer,  near  the  wall,  where  ambipolar  diffu¬ 
sion  is  significant,  and  an  outer  layer,  where  diffusion 
plays  no  significant  role.  The  first  subsection  following 
describes  the  model,  and  the  second  presents  results. 

5.1  Formulation 

We  start  with  G  =  — ^  =  constant  in  the  ionization 
zone.  The  differential  form  of  the  momentum  equation 
(7)  is  used: 


_dti  d 


p+ 


2n„ 


(35) 


where  the  pressure  for  a  two  temperature,  electrically 
neutral  plasma  is 


p  =  min,r|(a  +  (1  -  a)0)  (36) 

Here,  because  of  the  acceleration,  the  alternative 
version  for  the  ion  slip  flux  is  used,  so  that,  from  equa¬ 
tion  13,  we  can  define  the  flux,  T  =  — UeVi,  from 
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(37) 
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Tliis  equation,  along  with  the  following  three,  con¬ 
stitute  the  equations  of  motion  for  the  quasi-one- 
dimensional  isothermal  flow  problem.  Starting  with 
equation  35,  and  using  equations  36,  37,  the  momen¬ 
tum  equation  becomes  [9] 


(38) 


(39) 


dx  1  _  a(i  _  e)  -  ^ 

Equation  6  may  be  rewritten  as: 

dr  ^da  CaOL  ^  . 

dx  ~  ^dx  nr  ~  ~  ®) 

and  the  final  equation  is  the  magnetic  field  equation, 
equation  9. 

Nondimensionalizing  these  equations  will  allow  us 
to  analyze  them  parametricsJly.  Defining  reference 

^’^*5  =  ^  = 

the  new  variables  are  u  =  T  =  ^ 

E  =  E/{urtfho)  and  i  =  The  resulting  nondi- 
mensional  equations  of  motion  are: 
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d-t  _  da  ,  _  A(l-a)'\ 


d4 


E  —  uh 


1  +  (S-1)9 

where  for  convenience  we  define 


(43) 


/  = 


l  +  a-l)3 


-(i) 


-1 


E  ^  uh 

In  addition  to  the  parameters  defined  previously, 

{Ca)tlo<TretUrtf  _  (i<j)|| 


€D  = 


Afn 


A  = 


SaiG 


u;^^fio<rTef  La 


The  first  far  downstream  (^  — ►  c»)  boundary  con¬ 
dition  will  be  that  E  =  fib,  which  implies  an  infinite 
length  and  hence  an  infinite  magnetic  Reynolds  num¬ 
ber.  This  motivates  us  to  choose  the  magnetic  field 
itself  as  the  independent  variable.  Dividing  the  three 


other  equations  of  motion  (40,  41,  42)  by  equation  43, 
we  arrive  at  the  following  set  of  equations: 

di  («) 

The  second  downstream  (00)  boundary  condition  is 
the  same  balance  between  transverse  ambipolar  diffu¬ 
sion  and  ionization  used  in  the  constant  speed  case, 
which  is  now: 


1  “  QJoo 


Cgul^  _  €,,5^ 


(47) 


Additional  boundary  conditions  are  the  conditions 
at  the  injector  wall,  and  the  internal  boundary  condi¬ 
tions  at  the  sonic  passage  point.  The  injector  wall  con¬ 
dition  (25)  still  holds,  except  in  the  nondimensional 
form  it  is  now 


7«  =  a.  (1  +  ^)  (48) 

where  the  speed  at  the  inlet  wall,  tio,  is  now  a  result 
of  the  calculation,  rather  than  a  freely  chosen  value, 
as  was  the  case  in  the  constant  speed  problem. 

If  subsonic  injection  is  assumed,  then  smooth  pas¬ 
sage  through  a  “sonic  point”  [3]  is  required  in  the 
steady  state.  The  sonic  point  is  characterized  as  a 
singularity  that  occurs  when  the  denominator  of  — 
(equation  44)  is  identically  zero  at  some  b  =  b,.  In 
order  to  be  physically  possible,  the  numerator  must 
be  zero  at  the  same  This  poses  a  difficulty  in  the 
full  set  of  equations,  since  trajectories  diverge  rapidly 
near  this  singularity. 

Therefore  an  “outer”  set  of  equations  is  used  “far” 
from  the  wall  -  and  the  sonic  point  is  found  to  be 
embedded  in  this  outer  layer.  Since  1,  diffu¬ 

sion  is  dropped  (^  «  0)  and  the  combination  7/ei>  is 
eliminated  from  equations  44  and  45,  arriving  at  the 
nondiffusive  layer  equations: 


*  15  (1+21^) 

(“) 

and  the  ion  slip  flux  is,  fxom  post-processing: 
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The  variables  at  sonic  point  (where  the  denomina¬ 
tor  of  equation  49  is  zero)  may  all  be  found  explicitly 
by  setting  both  the  numerator  and  denominator  of 
equation  49  equal  to  zero,  and  using  the  integrated 
momentum  equation  (eqn.  7,  evaluating  the  constant, 
F  once  E  and  too  are  chosen)  to  solve  for  ,  a, ,  and 
Ug,  These  non-diifusive  equations  are  used  to  march 
back  toward  the  wall  from  the  sonic  point  (and  also 
from  the  sonic  point  out  to  oo  for  the  magnetic  diffu¬ 
sive  layer)  until  the  slip  7  becomes  positive,  and  then 
patched  with  the  diffusive  solution,  equations  44  to 
46,  are  to  get  to  the  wall  itself.  The  diffusive  set  is 
required  to  meet  the  boundary  condition  on  the  slip. 

At  this  point,  we  have  solved  the  problem  from  the 
wall  {b  =  1)  to  the  end  of  the  magnetic  diffusion  layer 
(where  E  =  iioofroo)-  This  process  yields  a  family  of 
boo{E)  for  each  E^  each  of  which  satisfies  the  bound¬ 
ary  conditions.  Thus,  there  is  one  remaining  degree 
of  freedom.  The  problem  is  closed  via  an  idealized 
channel  downstream  (the  “outer”  magnetic  problem), 
where  the  flow  is  assumed  to  be  constant  temperature, 
frozen  (at  a  =  aoo)i  and  zero  current  {E  =  ub).  The 
magnetic  field  is  then  a  function  of  the  area  of  the 
channel,  and  we  look  for  the  throat.  [3]  Now  (switch¬ 
ing  to  dimensional  variables),  <l>  =  EH  is  a  constant, 
so  that  the  speed  and  density  {miUg  =  m/{uA))  are 


u  =  —=  — 

BE  ^  (i>  w 

and  the  pressure,  from  equation  36  is 
<p  w 

so  that  the  momentum  equation,  becomes 


_j - I  dp  4-  d- — 

2  miUg  \  2po 


We  then  find  the  throat,  given  boo  aud  Ey  by  eval¬ 
uating  this  expression  for  decreasing  b  until  Ho/H  is 
minimized.  In  practice,  the  contraction  ratio,  Ho/Ht 
is  taken  as  a  given,  and  equation  53  gives  another  re¬ 
lationship  between  E  and  600. 

This  closes  the  problem  for  a  given  set  of  the  param¬ 
eters.  The  numerical  integrations  are  solved  using  a 
Runge-Kutta  space-marching  scheme,  starting  at  the 
sonic  point  (non-diffusive  equations),  and  marching 
out  both  towards  600  and  The  diffusive  (equations 
44,  45,  and  46)  and  non-diffusive  (equations  49  and  50) 
solutions  are  patched  just  after  the  outer,  non-diffusive 
set  value  for  7  is  small  and  positive  (>  0.001),  and 
then  the  diffusive  inner  set  is  integrated  to  the  wall. 

5.2  Results 

Following  are  results  from  the  inner-outer  approach  to 
the  accelerating  ignition  problem  as  described  above. 
The  plots  vs  (  have  been  truncated  so  as  to  show 
the  ionization  layer  more  clearly.  A  “standard”  set 
of  parameters  is  defined  in  this  work  corresponding 
to  roughly  Tg  ^  20000iir,  H  =  0.02m,  G  «  7.5  x 
10^^  m-2/5:  A  =  2,  ei?  =  0.005,  A3  =  0.1,  0  =  0.05, 
6  =  0.10,  and  q  =  0.001.  Since  the  ionization  rate 
coefficient  varies  the  most  dramatically  with  temper¬ 
ature  of  the  parameters,  as  can  be  seen  in  figure  4,  A 
will  be  varied  while  the  other  parameters  are  fixed. 
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Figure  4:  The  parameter  A  vs.  T*. 


Using  the  expressions  for  a,  mingy  and  p  above,  and 
the  nondimensional  variables  and  parameters  already 
defined,  plus  the  nondimensional  potential, 

^  ErtfSt 


Figure  5  shows  7  vs  a,  for  the  standard  set  of  pa¬ 
rameters,  and  a  contraction  ratio  of  one  {E  =  0.6019, 
boo  =  0.6583),  The  diffusive  and  non-diffusive  solu¬ 
tions  are  patched  at  7  =  0.0036,  a  =  0.5592,  and  the 
sonic  point  is  at  5  =  0.9854,  a,  =  0.9029.  The  iso¬ 
clines  are  for  the  diffusive  inner  set  at  6  =  1,  using 


then  equation  52  may  be  recast  and  integrated  to  ob-  equation  7  to  solve  for  u(cr).  Note  that,  although  the 
tain  the  following  expression:  inner  equations  diverge  as  7  0  (as  can  be  seen  from 

the  iso<^es),  the  outer  set  passes  through  smoothly, 
and  the  two  sets  patch  quite  well  at  a  small  positive 
7,  as  can  also  be  seen  from  the  isoclines. 

Figure  6  shows  the  ionization  fraction,  a,  vs  the 
nondimensional  length,  ^  =  as/Am,  for  the  “standard” 


Ho  .  /  1  .  2/3(0 +  aoc(l-^))ln(V)  ,  4(6oo  -  *») 

(53) 
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